edapy_06

In [1]:

    # edapy_06_00 clear all variables and import various modules
    %reset -f
    import os
    from datetime import date
    from math import exp, pi, sin, sqrt, floor, ceil
    import numpy as np
    import scipy.linalg as la
    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;

In [2]:

    # edapy_06_01
    # cosine wave

    # make the auxillary variable t
    N=100;
    Dt=1;
    tmin=0.0;
    tmax = tmin+Dt*(N-1);
    t = eda_cvec( np.linspace(tmin,tmax,N) );

    # a cosine of ampmitude C, time offset t0 and period T
    t0 = 10.0;
    T = 50;
    C= 2.0;
    d = C*np.cos( 2*pi*(t-t0)/T );

    # plot data
    fig1 = plt.figure(1,figsize=(8,4));
    ax1 = plt.subplot(1,1,1);
    plt.axis([tmin, tmax, 1.1*np.min(d), 1.1*np.max(d)])
    plt.plot(t,d,'k-');
    plt.xlabel('time t');
    plt.ylabel('d(t)');
    plt.title('Delayed cosine wave');
    plt.show();

In [57]:

    # edapy_06_02
    # plot columns of G
    # (not normalized - see eda06_13 for properly normalized version)

    # square data kernel
    N=32; # must be an even number
    M=N;

    # make the auxillary variable t
    N=100;
    Dt = 1;
    tmin=0.0;
    tmax = tmin+Dt*(N-1);
    Dt = (tmax - tmin)/(N-1);
    t = eda_cvec( np.linspace(tmin,tmax,N) );

    # standard definition of nyquist frequency fmax
    # frequency increment Df, frequency vector f,
    # angular frequency vector w
    Nf=floor(N/2+1);
    fmax = 1/(2*Dt);
    Df = fmax/(N/2);
    f = eda_cvec( np.linspace(0,fmax,Nf) );
    Nw=Nf;
    wmax = 2*pi*fmax;
    Dw = wmax/(N/2);
    w = eda_cvec( np.linspace(0,wmax,Nw) );

    # set up data kernel G
    G = np.zeros((N,M));

    # zero frequency column
    G[:,0] = np.ones((N,1)).ravel();

    # interior M/2-1 columns.  Note indice arithmetic
    # to keep track of pairs of columns, one for cos(),
    # one for sin()
    Mo2 = floor(M/2);
    for i in range(1,Mo2):
        j = 2*i-1;
        k = j+1;
        G[:,j]=np.cos(w[i,0]*t).ravel();
        G[:,k]=np.sin(w[i,0]*t).ravel();

    # nyquist column
    G[:,M-1] = np.cos(w[Nw-1,0]*t).ravel();

    # plot selected columns of dat kernel G
    fig1 = plt.figure(1,figsize=(8,4));
    plt.rcParams['xtick.bottom'] = False;
    plt.rcParams['xtick.labelbottom'] = False;
    plt.rcParams['xtick.top'] = True;
    plt.rcParams['xtick.labeltop'] = True;
    ip = 1;
    for i in [0,1,2,3,4,M-1]:
        ax = plt.subplot(1,6,ip);
        ax.invert_yaxis();
        plt.axis([-2, 2, tmin, tmax ]);
        plt.plot(G[:,i],t,'k-');
        plt.plot(G[:,i],t,'k.');
        mystr = "col %d" % (i+1);
        plt.title(mystr);
        if (i==0):
            plt.ylabel('time t');
        ip=ip+1;
    plt.show();

In [58]:

    # edapy_06_03: example of aliasing

    # make the auxillary variable t
    N=100;
    Dt=1;
    tmin=0.0;
    tmax = tmin+Dt*(N-1);
    Dt1 = (tmax - tmin)/(N-1);
    t = eda_cvec( np.linspace(tmin,tmax,N) );

    # standard set-up of Fourier parameters
    # note: N required to be an even integer
    Nf=floor(N/2+1);  # number of non-negative frequencies
    fmax = 1/(2*Dt1); # Nyquist frequency
    Df = fmax/(N/2);  # frequency increment
    # non-negative frequency axis
    f = eda_cvec( np.linspace(0,fmax,Nf) );
    Nw=Nf;            # number of angular frequencie
    wmax = 2*pi*fmax; # Nyquist angukar frequency
    Dw = wmax/(N/2);  # angular frequency incremnt
    # non-negative angular frequency axis
    w = eda_cvec( np.linspace(0,wmax,Nw) );

    # low frequency oscillation
    w1= 2*Dw;
    print('low frequency', w1);
    d1 = np.cos( w1*t );

    # high frequency oscillation
    w2= (2+N)*Dw;
    print('high frequency', w2);
    d2 = np.cos( w2*t );

    # set up higher sample rate time
    N=500;
    Dt2=0.1;
    tmin=0.0;
    tmax = tmin+Dt2*(N-1);
    t3 = eda_cvec( np.linspace(tmin,tmax,N) );
    d3 = np.cos( w2*t3 );

    # plot signals
    fig1 = plt.figure(1,figsize=(16,4));
    plt.rcParams['xtick.bottom'] = False;
    plt.rcParams['xtick.labelbottom'] = False;
    plt.rcParams['xtick.top'] = True;
    plt.rcParams['xtick.labeltop'] = True;

    ax1 = plt.subplot(2,1,1);
    plt.axis([tmin, tmax, -2, 2 ]);
    plt.plot(t,d1,'k-');
    plt.plot(t,d2,'ko');
    plt.xlabel('time t');
    plt.ylabel('d(t)');

    ax1 = plt.subplot(2,1,2);
    plt.axis([tmin, tmax, -2, 2 ]);
    plt.plot(t,d2,'k-');
    plt.plot(t,d2,'ko');
    plt.plot(t3,d3,'k:');
    plt.xlabel('time t');
    plt.ylabel('d(t)');
    plt.show();

    low frequency 0.12566370614359174
    high frequency 6.408849013323178

In [2]:

    # edapy_06_04: another example of aliasing

    # make the auxillary variable t
    N=50;
    Dt=1;
    tmin=0.0;
    tmax = tmin+Dt*(N-1);
    Dt1 = (tmax - tmin)/(N-1);
    t = eda_cvec( np.linspace(tmin,tmax,N) );

    # standard set-up of Fourier parameters
    # note: N required to be an even integer
    Nf=floor(N/2+1);  # number of non-negative frequencies
    fmax = 1/(2*Dt1); # Nyquist frequency
    Df = fmax/(N/2);  # frequency increment
    # non-negative frequency axis
    f = eda_cvec( np.linspace(0,fmax,Nf) );
    Nw=Nf;            # number of angular frequencie
    wmax = 2*pi*fmax; # Nyquist angukar frequency
    Dw = wmax/(N/2);  # angular frequency incremnt
    # non-negative angular frequency axis
    w = eda_cvec( np.linspace(0,wmax,Nw) );

    # low frequency oscillation
    w1= -2*Dw;
    print('low frequency', w1);
    d1 = np.sin( w1*t );

    # high frequency oscillation
    w2= (N-2)*Dw;
    print('high frequency', w2);
    d2 = np.sin( w2*t );

    # set up higher sample rate time
    N=500;
    Dt2=0.1;
    tmin=0.0;
    tmax = tmin+Dt2*(N-1);
    t3 = eda_cvec( np.linspace(tmin,tmax,N) );
    d3 =  np.cos( w2*t3 );

    # plot signals
    fig1 = plt.figure(1,figsize=(8,4));
    plt.rcParams['xtick.bottom'] = False;
    plt.rcParams['xtick.labelbottom'] = False;
    plt.rcParams['xtick.top'] = True;
    plt.rcParams['xtick.labeltop'] = True;

    ax1 = plt.subplot(2,1,1);
    plt.axis([tmin, tmax, -2, 2 ]);
    plt.plot(t,d1,'k-');
    plt.plot(t,d2,'ko');
    plt.xlabel('time t');
    plt.ylabel('d(t)');

    ax1 = plt.subplot(2,1,2);
    plt.axis([tmin, tmax, -2, 2 ]);
    plt.plot(t,d2,'k-');
    plt.plot(t,d2,'ko');
    plt.plot(t3,d3,'k:');
    plt.xlabel('time t');
    plt.ylabel('d(t)');
    plt.show();

    low frequency -0.25132741228718347
    high frequency 6.031857894892403

In [8]:

    # edapy_06_05: a.s.d. of the Neuse River hydrograph
    # using least-squares
    # (not normalized - see edapy_06_14 for properly normalized version)

    # load data from file
    D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
    [Nraw, K]=D.shape;
    traw = eda_cvec( D[:,0] );  # traw for "t raw"
    Dt = traw[1,0]-traw[0,0];  # sampling interval
    draw = eda_cvec( D[:,1] ); # draw for "d raw"

    # round off to even number of points
    No2 = floor(Nraw/2);
    N = 2*No2;
    d = eda_cvec( draw[0:N,0] );
    t = eda_cvec( traw[0:N,0] );

    # plot data
    fig1 = plt.figure(1,figsize=(10,4));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(1,1,1);
    plt.axis([np.min(t), np.max(t), 0, 1.1*np.max(d) ]);
    plt.plot(t,d,'k-');
    plt.xlabel('time (days)');
    plt.ylabel('discharge (cfs)');
    plt.title('Neuse River Hydrograph');
    plt.show();

    # nyquist frequencies
    Nf=floor(N/2+1);
    fmax = 1/(2*Dt);
    Df = fmax/(N/2);
    f = np.zeros((Nf,1));
    f[:,0] = np.linspace(0,fmax,Nf);
    Nw=Nf;
    wmax = 2*pi*fmax;
    Dw = wmax/(N/2);
    w = np.zeros((Nw,1));
    w[:,0] = np.linspace(0,wmax,Nw);
    fpos=np.zeros((Nf,1));
    fpos[:,0] = Df * np.linspace(0,No2,Nf);
    wpos=np.zeros((Nf,1));
    wpos[:,0] = Dw * np.linspace(0,No2,Nf);

    # set up data kernel G
    M=N;
    Mo2=floor(M/2);
    G = np.zeros((N,M));
    G[0:N,0:1] = np.ones((N,1)); # zero frequency column
    # interior M/2-1 columns.
    for i in range(1,Mo2):
        j = 2*i-1;
        k = j+1;
        G[0:N,j:j+1]=np.cos(w[i,0]*t);
        G[0:N,k:k+1]=np.sin(w[i,0]*t);
    G[0:N,M-1:M] = np.cos(w[Nw-1,0]*t); # nyquist column

    # solve least squares problem using
    # analytic formula for inv(GT*G)
    GTd = np.matmul( G.T, d );
    gtgi = 2*np.ones((M,1))/N;
    gtgi[0,0]=1/N;
    gtgi[M-1,0]=1/N;
    mest = np.multiply( gtgi, GTd );

    # compute power spectral density s
    s=np.zeros((Nw,1));
    s[0,0]=np.sqrt(mest[0,0]**2); # zero frequency
    for i in range(1,Mo2): # interior points
        j = 2*i-1;
        k = j+1;
        s[i,0] = sqrt(mest[j,0]**2 + mest[k,0]**2);
    s[Nw-1,0]= sqrt(mest[M-1,0]**2); # Nyquist frequency

    # plot spectrum
    fig2 = plt.figure(2,figsize=(10,4));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(1,3,1);
    plt.axis([0, fmax, 0, 1.1*np.max(s[1:Nf,0]) ]);
    plt.plot(fpos[1:Nf,0],s[1:Nf,0],'k-');
    plt.xlabel('frequency f (cycles per day)');
    plt.ylabel('amplitude spectrum');
    plt.title('linear in frequency');

    ax1 = plt.subplot(1,3,2);
    Tmax = 1000;
    plt.axis([0, Tmax, 0, 1.1*np.max(s[1:Nf,0]) ]);
    plt.plot(np.reciprocal(fpos[1:Nf,0]),s[1:Nf,0],'k-');
    plt.xlabel('periot T (days)');
    plt.ylabel('amplitude spectrum');
    plt.title('linear in period');

    ax1 = plt.subplot(1,3,3);
    plt.semilogx(fpos[1:Nf,0], s[1:Nf,0], 'k-', basex=10);
    plt.axis([0.001, 1, 0, 1.1*np.max(s[1:Nf,0])])
    plt.xlabel('frequency f (cycles per day)');
    plt.ylabel('amplitude spectrum');
    plt.title('logarithmic in frequency');

    plt.show();

    # check results
    dpre =np.matmul( G, mest );
    e = d - dpre; # individual errors
    E = np.matmul( e.T, e ); # total error
    print('Error', E);

    Error [[2.87461543e-15]]

In [14]:

    # edapy_06_06: a.s.d. of the Neuse River hydrograph
    # using fft command
    # (unnormalized spectral density - see eda06_14 for proper normalization)

    # load data from file
    D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
    [Nraw, K]=D.shape;
    traw = eda_cvec( D[:,0] ); # traw for "t raw"
    Dt = traw[1,0]-traw[0,0]; # sampling interval
    draw = eda_cvec( D[:,1] ); # draw fro "d raw"

    # round off to even number of points
    No2 = floor(Nraw/2);
    N=2*No2;
    d = eda_cvec( draw[0:N,0] ); # truncate data vector
    t = eda_cvec( traw[0:N,0] ); # truncate time vector

    # plot data
    fig1 = plt.figure(1,figsize=(10,4));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(1,1,1);
    plt.axis([np.min(t), np.max(t), 0, 1.1*np.max(d) ]);
    plt.plot(t,d,'k-');
    plt.xlabel('time (days)');
    plt.ylabel('discharge (cfs)');
    plt.title('Neuse River Hydrograph');
    plt.show();

    No2 = floor(N/2); # N/2 rounded off to an integer
    M=N; # number of model parameters equals number of data

    # set up time axis
    tmin = 0;
    tmax = Dt*(N-1);
    t= eda_cvec( np.linspace(tmin,tmax,N) );

    # generic frequency setup
    fmax=1/(2.0*Dt);  # nyquist frequency
    Df=fmax/No2;      # frequency sampling
    Nf=No2+1;         # number of non-negative frequencies
    # f is full freqeuncy vector:
    # f = [0, Df, 2Df ... fmax, -fmax+Df, ... -Df]
    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 angular frequency vector
    Nw=Nf;            # number of non-negative angular f's
    fpos = eda_cvec(Df * np.linspace(0,No2,Nf)); # non-neg f's
    wpos=2*pi*fpos;   # non-negative angular frequencies

    # Compute fourier coefficients. Be careful not to forget
    # the axis=0 (fft over elements of column). The default
    # axis=1 (fft over elements of rows) won't work!
    mest = eda_cvec( np.fft.fft(d, axis=0) );

    # ampltude spectrum (note is un-normalized)
    mestpos = eda_cvec( mest[0:Nf,0] ); # pos f's only
    s = np.abs(mestpos);                # a.s.d.

    # plot spectrum
    fig2 = plt.figure(2,figsize=(10,4));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    # linear in frequency plot
    ax1 = plt.subplot(1,3,1);
    plt.axis([0, fmax, 0, 1.1*np.max(s[1:Nf,0]) ]);
    plt.plot(fpos[1:Nf,0],s[1:Nf,0],'k-');
    plt.xlabel('frequency f (cycles per day)');
    plt.ylabel('amplitude spectrum');
    plt.title('linear in frequency');

    # linear in perior plot
    ax1 = plt.subplot(1,3,2);
    Tmax = 1000;
    plt.axis([0, Tmax, 0, 1.1*np.max(s[1:Nf,0]) ]);
    plt.plot(np.reciprocal(fpos[1:Nf,0]),s[1:Nf,0],'k-');
    plt.xlabel('periot T (days)');
    plt.ylabel('amplitude spectrum');
    plt.title('linear in period');

    # semilog plot
    ax1 = plt.subplot(1,3,3);
    plt.semilogx(fpos[1:Nf,0], s[1:Nf,0], 'k-', basex=10);
    plt.axis([0.001, 1, 0, 1.1*np.max(s[1:Nf,0])])
    plt.xlabel('frequency f (cycles per day)');
    plt.ylabel('amplitude spectrum');
    plt.title('logarithmic in frequency');

    plt.show();

    # compute Error of d->mest->dnew
    # Inverse FFT: Don't forget the axis=0!
    dnew = eda_cvec( np.real(np.fft.ifft(mest,axis=0)));
    e = d - dnew; # individual errors
    E = np.matmul( e.T, e );  # total error
    print('Error of d->mest->dnew:', E );

    Error of d->mest->dnew: [[1.92672359e-20]]

In [62]:

    # edapy_06_07
    # spectra of Normal curves of varying widths

    # time series of length N=256 with sampling interval of Dt=1

    # generic time set up
    N=256;   # must be an even number
    Dt=1.0;
    tmin = 0;
    tmax = Dt*(N-1);
    t = eda_cvec( np.linspace(tmin,tmax,N) );

    # generic frequency set up
    No2 = floor(N/2);  # N divided by two as an integer
    fmax=1/(2.0*Dt);   # nyquist frequency
    Df=fmax/No2;       # frequency increment
    Nf=No2+1;          # number of non-negative frequencies
    # full vector of frequencies
    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 vector of angular fequencies
    Nw=Nf;             # number of angular frequencies
    fpos=eda_cvec(Df*np.linspace(0,No2,Nf)); # non-negative frequencies
    wpos=2*pi*fpos;    # non-negative angular frequencies

    # stamdard deviations, chosen to make an aestheic plot
    Nsd=6;
    sd=eda_cvec( [1.25, 2.5, 5, 10, 20, 40] );

    # matrix D with each column a time series, one column for each sd
    D = np.zeros((N,Nsd));
    # matrix Ds with each column an amplitude spectral density
    Ds = np.zeros((Nf,Nsd));

    # populate D and Ds
    for j in range(Nsd):
        t1 = tmax/2;     # mean
        sd2 = sd[j]**2;  # variance
        d = np.exp( -np.power(t-t1,2) / (2*sd2) );  # normal function
        D[:,j] = d.ravel();

        #  compute fourier coefficients
        dfft = eda_cvec( np.fft.fft(d, axis=0) );

        # delete negative frequencies
        dfftpos = eda_cvec( dfft[0:Nf,0] );

        # ampltude spectrum
        Ds[:,j]=np.abs(dfftpos).ravel();

    eda_draw(D[:,0],D[:,1],D[:,2],D[:,3],D[:,4],D[:,5],'    ',Ds[:,0],Ds[:,1],Ds[:,2],Ds[:,3],Ds[:,4],Ds[:,5]);

In [63]:

    # edapy_06_08 spike function

    # generic time set up
    N=256;  # length of time sereis; must be even
    Dt=1.0; # samppling interval
    tmin = 0;
    tmax = Dt*(N-1);
    t = eda_cvec( np.linspace(tmin,tmax,N) ); # time vector

    # spike time series of length N=256 with sampling interval of Dt=1
    d = np.zeros((N,1));
    d[0,0]=1.0;

    # generic frequency setup
    No2 = floor(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 frequencies
    Nw=Nf; # number rof non-negative angular frequencies
    fpos = eda_cvec(Df * np.linspace(0,No2,Nf) ); # non-negative frequencies
    wpos= 2*pi*fpos; # non-negathive angular frequencies

    # take fourier transform
    # Note: Don't forget axis=0 option!
    dt = eda_cvec( np.fft.fft(d, axis=0) );

    # plot
    fig1 = plt.figure(1,figsize=(10,4));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(2,1,1);
    plt.axis( [-tmax/10, tmax, -0.5, 1.5 ]);
    plt.plot(t,d,'k-');
    plt.plot(t,d,'k.');
    plt.xlabel('time t');
    plt.ylabel('d(t)');
    plt.title('time series');

    ax1 = plt.subplot(2,1,2);
    plt.axis( [-fmax, fmax, 0, 2] );
    plt.plot(f,np.real(dt),'k-');  # the transform should be purely real, so take real part
    plt.plot(f,np.real(dt),'k.');
    plt.xlabel('frequency f');
    plt.ylabel('d-tilde(f)');
    plt.title('time series');

    plt.show();

In [64]:

    # edapy_06_09: area under curve using FFT

    # time series of length N=256 with sampling interval of Dt=1

    # generic timeset up
    N=256; # length of time series, must be even
    Dt=1.0; # sampling interval
    tmin = 0; # minimum time
    tmax = Dt*(N-1); # maximum time
    t = eda_cvec( np.linspace(tmin,tmax,N) ); # time vector

    # generic frequency up
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # nyquist frequency
    Df=fmax/No2;       # frequency smaping
    Nf=No2+1;          # number on non-negative frequencies
    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 frequeny increment
    w=2*pi*f;          # full (pos and neg) ana
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies
    wpos=2*pi*fpos;    # non-negative angular frequencies

    # exemplary function
    w0 = 2*pi*fmax/20;  # frequency of oscilation
    d= np.multiply(np.cos(w0*t),np.exp(-0.25*w0*t)); # product of cos(0 and delining exp()

    # Fourier transform; note normalization factor of Dt
    # and use of axis=0 option.
    dt = eda_cvec( Dt*np.fft.fft(d, axis=0) );
    # area by fft
    area2 = np.real(dt[0,0]);

    # area by Reimann summation
    area1 = Dt*np.sum(d);

    # plot
    fig1 = plt.figure(1,figsize=(10,4));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(1,1,1);
    plt.axis( [0, tmax, -1, 1] );
    plt.plot(t,d,'k-');
    plt.xlabel('time t');
    plt.ylabel('d(t)');
    titlestr = "sum area %f     fft area %f" % (area1, area2);
    plt.title(titlestr);

    plt.show();

In [16]:

    # edapy_06_10: time shift using FFT

    # generic time= set up
    N=256;
    Dt=1.0;
    tmin = 0.0;
    tmax = Dt*(N-1);
    t=eda_cvec(np.linspace(tmin,tmax,N));

    # generic time/frequency set up
    No2 = floor(N/2);
    fmax=1/(2.0*Dt);
    Df=fmax/No2;
    Nf=No2+1;
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;
    w=2*pi*f;
    Nw=Nf;
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf));
    wpos=2*pi*fpos;

    # exemplary function
    d = np.zeros((N,1));
    w0 = 2*pi*fmax/20;
    d[:,0] = np.multiply(np.cos(w0*t),np.exp(-0.25*w0*t)).ravel();

    # time shift
    t0 = t[15,0];
    dt = (np.multiply(
        np.exp((complex(0,-1)*w*t0)), np.fft.fft(d, axis=0)));

    # inverse Fourier transform
    # result should be purely real, so take real part
    # to get rid of zero imaginary part
    ds = eda_cvec( np.real((np.fft.ifft(dt,axis=0))));

    # plot
    fig1 = plt.figure(1,figsize=(10,7));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(2,1,1);
    plt.axis( [0, tmax, -1, 1] );
    plt.plot(t,d,'k-');
    plt.ylabel('d(t)');
    plt.title('time series before shift');

    ax1 = plt.subplot(2,1,2);
    plt.axis( [0, tmax, -1, 1] );
    plt.plot(t,ds,'k-');
    plt.xlabel('time t');
    plt.ylabel('d(t)');
    titlestr = "time series after rshift of %.2f" % (t0);
    plt.title(titlestr);

    plt.show();

    (256, 1)

In [19]:

    # eda06_11: derivative using FFT

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

    # generic time/frequency set up
    No2 = floor(N/2);
    fmax=1/(2.0*Dt);
    Df=fmax/No2;
    Nf=No2+1;
    f=eda_cvec(Df * np.concatenate((np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));
    Dw=2*pi*Df;
    w=2*pi*f;
    Nw=Nf;
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf));
    wpos=2*pi*wpos;

    # exemplary function
    d = np.zeros((N,1));
    t0 = tmax/2;
    sd = 50;
    d[:,0] = np.exp(-np.power((t-t0),2)/(2*sd**2)).ravel();

    # take Fourier transform and multiply by iw
    # except for zero freuency
    dt = np.multiply(complex(0,1)*w , np.fft.fft(d, axis=0));
    # inverse Fourier transform
    dddt = np.real((np.fft.ifft(dt,axis=0)));

    # derivative by finite differences
    dddt2 = np.zeros((N,1));
    dddt2[0:N-1,0]=(d[1:N,0]-d[0:N-1,0])/Dt;

    # plot
    fig1 = plt.figure(1,figsize=(10,8));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(3,1,1);
    plt.axis( [0, tmax, -1.1*np.max(d), 1.1*np.max(d)] );
    plt.plot(t,d,'k-');
    plt.ylabel('d(t)');

    ax1 = plt.subplot(3,1,2);
    plt.axis( [0, tmax, -1.1*np.max(dddt), 1.1*np.max(dddt)] );
    plt.plot(t,dddt,'k-');
    plt.xlabel('time t');
    plt.ylabel('dd/dt by fft');

    ax1 = plt.subplot(3,1,3);
    plt.axis( [0, tmax, -1.1*np.max(dddt2), 1.1*np.max(dddt2)] );
    plt.plot(t,dddt2,'k-');
    plt.xlabel('time t');
    plt.ylabel('dd/dt by sum');
    plt.show();

In [24]:

    # edapy_06_12: integral using FFT

    # generic time/y set up
    N=256;
    Dt=1.0;
    tmin = 0;
    tmax = Dt*(N-1);
    t=eda_cvec(np.linspace(tmin,tmax,N));

    # generic frequency set up
    No2 = floor(N/2);
    fmax=1/(2.0*Dt);
    Df=fmax/No2;
    Nf=No2+1;
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;
    w=2*pi*f;
    Nw=Nf;
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf));
    wpos=2*pi*fpos;

    # exemplary function
    sd = 50;
    t1 = tmax/4;
    t2 = 3*tmax/4;
    sd = 25;
    d = (np.exp(-np.power((t-t1),2)/(2*sd**2)) - np.exp(-np.power((t-t2),2)/(2*sd**2))  );

    # take Fourier transform and divide by iw
    # except set zero frequency to zero
    wr=np.concatenate(
     (np.zeros((1,1)),np.reciprocal(complex(0,1)*w[1:N,0:1])) );
    dt = np.multiply(wr, np.fft.fft(d, axis=0));

    # inverse Fourier transform, don't forget the axis=0
    # since result should be purely real, delete imaginary part
    integral1 = eda_cvec( np.real((np.fft.ifft(dt,axis=0))));

    # integeral by Reimann sum
    integral2 = eda_cvec((Dt*np.cumsum(d)));

    # adjust so that integrals have the same integration constant (zero)
    integral1 = integral1 - integral1[0,0];
    integral2 = integral2 - integral2[0,0];

    # plot
    fig1 = plt.figure(1,figsize=(10,8));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(3,1,1);
    plt.axis( [0, tmax, -1.1*np.max(d), 1.1*np.max(d)] );
    plt.plot(t,d,'k-');
    plt.ylabel('d(t)');

    ax1 = plt.subplot(3,1,2);
    plt.axis( [0, tmax, -1.1*np.max(integral1), 1.1*np.max(integral1)] );
    plt.plot(t,integral1,'k-');
    plt.xlabel('time t');
    plt.ylabel('integral by fft');

    ax1 = plt.subplot(3,1,3);
    plt.axis( [0, tmax, -1.1*np.max(integral2), 1.1*np.max(integral2)] );
    plt.plot(t,integral2,'k-');
    plt.xlabel('time t');
    plt.ylabel('dd/dt by diff');
    plt.show();

In [25]:

    # edapy_06_13: Crib Sheet 6.2
    # makes exemplary time series d(t),
    # amplitude spectral density s(f) pairs for

    # time series of length N=256 with sampling interval of Dt=1
    N=128;
    Dt=0.5;
    tmin = 0;
    tmax = Dt*(N-1);
    t=eda_cvec(np.linspace(tmin,tmax,N));

    # generic frequency set up
    No2 = floor(N/2);
    fmax=1/(2.0*Dt);
    Df=fmax/No2;
    Nf=No2+1;
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;
    w=2*pi*f;
    Nw=Nf;
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf));
    wpos=2*pi*fpos;

    # duration of time series
    T = tmax-tmin;

    # plot 1
    fig1 = plt.figure(1,figsize=(10,8));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;
    p=1;

    # narrow function
    sigmad = 1.0;
    d = (np.exp(-np.power(t-tmax/2,2)/(2*(sigmad**2)))+np.random.normal(0,0.1,(N,1)));

    dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );   # fourier transform, note normalization of Dt
    dfftpos = eda_cvec( dfft[0:Nf,0] ) ; # positive fequencies, only
    ds=sqrt(2/T)*np.abs(dfftpos); # properly normalized amplitude spectral density

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, T, -1.1, 1.1] );
    plt.plot(t,d,'k-');
    plt.ylabel('d(t)');
    p=p+1;

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
    plt.plot(fpos,ds,'k-');
    plt.ylabel('s(f)');
    p=p+1;

    # zero area wavelet
    sigmad1 = 2.5;
    sigmad2 = 1.0;
    d = np.multiply( np.exp(-np.power(t-T/2,2)/(2*(sigmad1**2))) , np.cos((t-T/2)/sigmad2));
    d = d + np.random.normal(0,0.1,(N,1));

    dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );
    dfftpos = eda_cvec( dfft[0:Nf,0] );
    ds=sqrt(2/T)*np.abs(dfftpos);

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, T, -1.1, 1.1] );
    plt.plot(t,d,'k-');
    plt.ylabel('d(t)');
    p=p+1;

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
    plt.plot(fpos,ds,'k-');
    plt.ylabel('s(f)');
    p=p+1;

    # wide function
    sigmad = 7.0;
    d = np.exp(-np.power(t-tmax/2,2)/(2*(sigmad**2))) + np.random.normal(0,0.1,(N,1));

    # fourier transform
    dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );  # note Dt normalization
    dfftpos = eda_cvec( dfft[0:Nf,0] );  # keep non-negative frequencies
    ds = sqrt(2/T)*np.abs(dfftpos); # amplitude spectral density

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, T, -1.1, 1.1] );
    plt.plot(t,d,'k-');
    plt.ylabel('d(t)');
    p=p+1;

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
    plt.plot(fpos,ds,'k-');
    plt.ylabel('s(f)');
    p=p+1;

    # triangle
    No4 = floor(N/4);
    d = np.zeros((N,1));
    d[No4:2*No4,0] = np.linspace(0,1,No4);
    d[2*No4:3*No4,0] = np.linspace(1,0,No4);
    d[:,0] = (d + np.random.normal(0,0.01,(N,1))).ravel();

    dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform with Dt normalization
    dfftpos = eda_cvec( dfft[0:Nf,0] );  # keep positive frequencies
    ds = sqrt(2/T)*np.abs(dfftpos);  # amplitude spectral density

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, T, -1.1, 1.1] );
    plt.plot(t,d,'k-');
    plt.xlabel('time t');
    plt.ylabel('d(t)');
    p=p+1;

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
    plt.plot(fpos,ds,'k-');
    plt.xlabel('frequency f');
    plt.ylabel('s(f)');
    p=p+1;

    plt.show();

    print("\n\n\n\n");

    # plot 2
    fig2 = plt.figure(2,figsize=(10,8));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;
    p=1;

    # approximately sinusoidal, with time varying frequency w0
    w0 = eda_cvec( (0.25*2*pi*(0.95+0.1*np.linspace(0,N-1,N)/(N-1))) );
    d = np.sin(np.multiply(w0,t)) + np.random.normal(0,0.2,(N,1));

    dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );
    dfftpos = eda_cvec( dfft[0:Nf,0] );
    ds = sqrt(2/T)*np.abs(dfftpos);

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, T, -1.1, 1.1] );
    plt.plot(t,d,'k-');
    plt.ylabel('d(t)');
    p=p+1;

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
    plt.plot(fpos,ds,'k-');
    plt.ylabel('s(f)');
    p=p+1;

    # beating, two sinusoids with frequencies in 2:1 ratio
    w1 = 0.25*2*pi;
    w2 = w1/2;
    d = np.sin(w1*t) + 0.5*np.sin(w2*t) + np.random.normal(0,0.2,(N,1));

    dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier trasnform with Dt normailzation
    dfftpos = eda_cvec( dfft[0:Nf,0] ); # non-negative fequencies only
    ds =sqrt(2/T)*np.abs(dfftpos).ravel(); # amplitude spectral density

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, T, -1.1, 1.1] );
    plt.plot(t,d,'k-');
    plt.ylabel('d(t)');
    p=p+1;

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
    plt.plot(fpos,ds,'k-');
    plt.ylabel('s(f)');
    p=p+1;

    # chirp, frequency that increases with time
    w0 = eda_cvec(0.25*2*pi*(2*np.linspace(0,N-1,N)/(N-1)));
    d = np.sin(np.multiply(w0,t)) + np.random.normal(0,0.05,(N,1));

    dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform with Dt normaization
    dfftpos = eda_cvec( dfft[0:Nf,0] );  # non-negative frequencies only
    ds = sqrt(2/T)*np.abs(dfftpos);

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, T, -1.1, 1.1] );
    plt.plot(t,d,'k-');
    plt.ylabel('d(t)');
    p=p+1;

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
    plt.plot(fpos,ds,'k-');
    plt.ylabel('s(f)');
    p=p+1;

    # sawtooth
    d = ((np.mod(t,N/10)/(N/10)-0.5)+np.random.normal(0.0,0.005,(N,1)));

    dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform with Dt normalization
    dfftpos = eda_cvec( dfft[0:Nf,0] ); # non-negative frequencies
    ds = sqrt(2/T)*np.abs(dfftpos); # amplitude spectral ratio

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, T, -1.1, 1.1] );
    plt.plot(t,d,'k-');
    plt.xlabel('time t');
    plt.ylabel('d(t)');
    p=p+1;

    ax1 = plt.subplot(4,2,p);
    plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
    plt.plot(fpos,ds,'k-');
    plt.xlabel('frequency f');
    plt.ylabel('s(f)');
    p=p+1;

    plt.show();

In [26]:

    # edapy_06_14: power spectral density
    # comparison of variances between time series and its power spectral density

    # random time series with zero mean and variance sigmad^2
    N=256;   # length of time series, must be even
    Dt=0.01; # sampling interval
    tmin = 0;
    tmax = Dt*(N-1);
    t=eda_cvec(Dt*np.linspace(0,N-1,N));
    sigmad1=10; # standard deviation
    sigmad1sq=sigmad1**2; # variance
    d = np.random.normal(0,sigmad1,(N,1)); # random timeseroes
    # this code will also work for a timeseries with non-zero mean,
    # but then the "variance" is really "power"
    d = d-np.mean(d);  # this statement can be commented out for "power"

    # variance estimated from time series
    sigmad2sq = np.matmul(d.T,d)/N;

    # standard frequency setup for cos() and sin() representation of data
    Nf=floor(N/2+1); # number o fnon-negative frequencies
    fmax = 1/(2*Dt); # nyquist frequency
    Df = fmax/(N/2); # frequency interval
    f = eda_cvec( np.linspace(0,fmax,Nf) ); # positive frequencies
    Nw=Nf; # number of non-negative angular frequencies
    wmax = 2*pi*fmax; # maximum angular frequency
    Dw = 2*pi*Df; # angular frequency interval
    w = 2*pi*f; # non-negative angular frequencies

    # set up data kernel G for fourier coefficients
    M=N;
    G = np.zeros((N,M));

    # zero frequency column
    G[:,0] = np.ones((N,1)).ravel();

    # interior M/2-1 columns
    # note index arithmetic due to columns occuring in pairs
    Mo2 = floor(M/2);
    for i in range(1,Mo2):
        j = 2*i-1;
        k = j+1;
        G[:,j]=np.cos(w[i,0]*t).ravel();
        G[:,k]=np.sin(w[i,0]*t).ravel();

    # nyquist column
    G[:,M-1] = np.cos(w[Nw-1,0]*t).ravel();

    # solve least squares problem using
    # analytic formula for inv(GT*G)
    GTd = np.matmul( G.T, d );
    # since inv( GT G ) is known analytically, use shortcut
    gtgi = 2*np.ones((M,1))/N;
    gtgi[0,0]=1/N;
    gtgi[M-1,0]=1/N;
    mest = np.multiply( gtgi , GTd );

    e = d - np.matmul(G,mest); # individual errors
    E = np.matmul( e.T , e ); # overall erros
    # print('error', E);  # uncomment to check that error is zero

    # generic frequency setup for fourier transform
    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, neg) frequencies
    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 freuqncy samppling
    w=2*pi*f;       # angular frequencies
    Nw=Nf;          # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies
    wpos=2*pi*fpos; # non-negative angular frequencies

    # discrete fourier transform, nno normalization
    C = eda_cvec( np.fft.fft(d, axis=0) );

    # inverse Fourier transform; result should be purely real, so take real part
    drecon = eda_cvec( np.real((np.fft.ifft(C,axis=0))) );
    e = d - drecon; # individual errors
    E = np.matmul( e.T, e ); # total error
    # print('error', E); # uncomment to check that error is zero

    # comparison of coefficients
    jf=10;
    kf=2*jf-1;
    printstr = "check agreement between sines/cosine and fft at frequency %d" % (jf);
    print(printstr);
    printstr = "A2 B2                      %.4f %.4f" % (mest[kf,0], mest[kf+1,0]);
    print(printstr);
    printstr = "(2/N)real(C) -(2/N)imag(C) %.4f %.4f" % ((2/N)*np.real(C[jf,0]), -(2/N)*np.imag(C[jf,0]));
    print(printstr);

    printstr = "variance estimated from time series            %.4f" % (sigmad2sq);
    print(printstr);

    sigmad3sq = 0.5*np.matmul(mest.T,mest);
    # normally one doesn't bother to correct the first and last frequency, but ...
    sigmad3sq = sigmad3sq  + 0.5*(mest[0,0]**2 + mest[N-1,0]**2)
    printstr = "variance estimated from sine/cosines           %.4f" % (sigmad3sq);
    print(printstr);

    # variance estimated from fft coefficients
    C2 = np.zeros((Nf,1))
    C2[:,0] = np.real(np.multiply(C[0:Nf,0].conj(),C[0:Nf,0])).ravel();
    sigmad4sq = (2/(N**2)) * np.sum(C2);
    # normally one doesn't bother to correct the first and last frequency, but ...
    sigmad4sq = sigmad4sq - 0.5*(2/(N**2))*C2[0,0] - 0.5*(2/(N**2))*C2[Nf-1,0];
    printstr = "variance estimated from fft coefficients       %.4f" % (sigmad4sq);
    print(printstr);

    # approximation to integral transform
    dbar=Dt*C;

    # variance estimated from fourier transform
    dbar2 = eda_cvec( np.real(np.multiply(dbar[0:Nf,0].conj(),dbar[0:Nf,0])));
    sigmad4sq = (2/((Dt**2)*(N**2))) * np.sum( dbar2 );
    # normally one doesn't bother to correct the first and last frequency, but ...
    sigmad4sq = sigmad4sq - 0.5*(2/((Dt**2)*(N**2)))*dbar2[0,0] - 0.5*(2/((Dt**2)*(N**2)))*dbar2[Nf-1,0];
    printstr = "variance estimated from fourier transform      %.4f" % (sigmad4sq);
    print(printstr);

    # integral from 0 to fny of |dbar(f)|2
    theintegral = Df * np.sum( dbar2 );
    sigmad5sq = (2/(Df*(Dt**2)*(N**2))) * theintegral;
    # normally one doesn't bother to correct the first and last frequency, but ...
    sigmad5sq = sigmad5sq - 0.5*Df*(2/(Df*(Dt**2)*(N**2)))*dbar2[0,0] - 0.5*Df*(2/(Df*(Dt**2)*(N**2)))*dbar2[Nf-1,0];
    printstr = "variance estimated from frequency integral     %.4f" % (sigmad5sq);
    print(printstr);

    # definition of power spectral density
    T = tmax-tmin+Dt;  # note the correction of the length of the window by Dt
    thepsd = (2/T) * dbar2;

    # integral of power spectral density
    sigmad6sq = Df * np.sum(thepsd);
    # normally one doesn't bother to correct the first and last frequency, but ...
    sigmad6sq  = sigmad6sq  - 0.5*Df*thepsd[0,0] - 0.5*Df*thepsd[Nf-1,0];
    printstr = "variance estimated from power spectral density %.4f" % (sigmad6sq);
    print(printstr);

    check agreement between sines/cosine and fft at frequency 10
    A2 B2                      -1.3605 0.5448
    (2/N)real(C) -(2/N)imag(C) -1.3605 0.5448
    variance estimated from time series            112.9074
    variance estimated from sine/cosines           112.9074
    variance estimated from fft coefficients       112.9074
    variance estimated from fourier transform      112.9074
    variance estimated from frequency integral     112.9074
    variance estimated from power spectral density 112.9074

In [30]:

    # edapy_06_15: p.s.d. of the Neuse River hydrograph
    # properly normalized amplitude spectral density of
    # Neuse River hydrograph using FFT

    # load data from file
    D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
    [Nraw, K]=D.shape;
    traw = eda_cvec( D[:,0] ); # time
    Dt = traw[1,0]-traw[0,0];  # sampling interval
    draw = eda_cvec( D[:,1] ); # data

    # round off to even number of points
    No2 = floor(Nraw/2);
    N=2*No2;
    d = eda_cvec( draw[0:N,0] );
    t = eda_cvec( traw[0:N,0] );

    # generic time set up
    tmin = 0;
    tmax = Dt*(N-1);
    t=eda_cvec(np.linspace(tmin,tmax,N));

    # generic frequency set up
    No2 = floor(N/2);  # integer version of N/2
    M=N;               # number of model parameters equal number of data
    fmax=1/(2.0*Dt);   # nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequenices
    # full (pos, 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, neg) anguar frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies
    wpos=2*pi*fpos;    # non-negative angular frequencies

    #  compute fourier coefficients. Note normalization of Dt
    dbar = Dt*np.fft.fft(d, axis=0); # fourier transform
    dbarpos = dbar[0:Nf,0:1]; # non-negative frequencies
    T = N*Dt;  # time series length
    # p.s.d.
    s2 = (2/T) * np.real(np.multiply(dbarpos.conj(), dbarpos));

    # plot data
    fig1 = plt.figure(1,figsize=(10,4));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(1,1,1);
    plt.axis([np.min(t), np.max(t), 0, 1.1*np.max(d) ]);
    plt.plot(t,d,'k-');
    plt.xlabel('time (days)');
    plt.ylabel('discharge (cfs)');
    plt.title('Neuse River Hydrograph');
    plt.show();

    # plot power spectral density
    fig2 = plt.figure(2,figsize=(10,4));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(1,3,1);
    plt.axis([0, fmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);
    plt.plot(fpos[1:Nf,0],s2[1:Nf,0],'k-');
    plt.xlabel('frequency f (cycles per day)');
    plt.ylabel('power spectral density');
    plt.title('linear in frequency');

    ax1 = plt.subplot(1,3,2);
    Tmax = 1000;
    plt.axis([0, Tmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);
    plt.plot(np.reciprocal(fpos[1:Nf,0]),s2[1:Nf,0],'k-');
    plt.xlabel('periot T (days)');
    plt.ylabel('power spectral density');
    plt.title('linear in period');

    ax1 = plt.subplot(1,3,3);
    plt.axis([0.001, 1, 0, 1.1*np.max(s2[1:Nf,0])]);
    plt.semilogx(fpos[1:Nf,0], s2[1:Nf,0], 'k-', basex=10);
    plt.xlabel('frequency f (cycles per day)');
    plt.ylabel('power spectral density');
    plt.title('logarithmic in frequency');

    plt.show();

    # variance calculation
    v1 = np.std(d)**2;
    printstr = "variance estimated from time series            %.6e" % (v1);
    print(printstr);
    v2 = Df * np.sum(s2[1:Nf,0]);
    printstr = "variance estimated from power spectral density %.6e" % (v2);
    print(printstr);

    variance estimated from time series            8.396180e+06
    variance estimated from power spectral density 8.396180e+06

In [32]:

    # eda06_15
    # make a random noise dataset

    # Note that he noise is saved to mynoise.txt, not noise.txt
    # to avoid overwriting the standard distribution.  So change
    # the name if you want to overwrite noise.txt

    # generic time setup
    N=1024;          # length of time series, must be even
    Dt=1.0;          # sampling interval
    tmin = 0.0;      # start time
    tmax = Dt*(N-1); # end time
    t = eda_cvec(np.linspace(tmin,tmax,N)); # time vector

    # make vector d of random noise
    dbar=0.0;   # mean
    sigmad=1.0; # standard deviation
    d = np.random.normal( dbar, sigmad, (N,1) );

    # plot data
    fig1 = plt.figure(1,figsize=(10,4));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    ax1 = plt.subplot(1,1,1);
    r = dbar+3*sigmad;  # this range encloses 99% of data
    plt.axis([tmin, tmax, -r, r ]);
    plt.plot(t,d,'k-');
    plt.plot(t,d,'k.');
    plt.plot([tmin,tmax],[0,0],'k-');
    plt.xlabel('time');
    plt.ylabel('d(t)');
    titlestr = "Normally-distributed random timeseries witm mean %.2f and sigma %.2f" % (dbar, sigmad);
    plt.title(titlestr);
    plt.show();

    # save matrix to file of tab-separated values
    # the file name has been changed from "noise.txt" to "noise2.txt"
    # so as not to overwrite the existing file in the distributuon
    Dm = np.concatenate( (t,d),axis=1);
    np.savetxt("../Data/munoise.txt",Dm, delimiter="\t");

In [ ]: