How to use the mudpy.green.stdecimate function in mudpy

#    bandpass=[1./20,1./2]
                    #    fsample=1./Ess[0].stats.delta
                    #    Ess[ktrace].data=lfilt(Ess[ktrace].data,bandpass,fsample,2)
                    #    Nss[ktrace].data=lfilt(Nss[ktrace].data,bandpass,fsample,2)
                    #    Zss[ktrace].data=lfilt(Zss[ktrace].data,bandpass,fsample,2)
                    #    Eds[ktrace].data=lfilt(Eds[ktrace].data,bandpass,fsample,2)
                    #    Nds[ktrace].data=lfilt(Nds[ktrace].data,bandpass,fsample,2)
                    #    Zds[ktrace].data=lfilt(Zds[ktrace].data,bandpass,fsample,2)
                    #    #bandpass=None
                    ### END AWFUL TERRIBLE REALLY VERY BAD HACK
                        
                    if decimate!=None: 
                        Ess[ktrace]=stdecimate(Ess[ktrace],decimate)
                        Nss[ktrace]=stdecimate(Nss[ktrace],decimate)
                        Zss[ktrace]=stdecimate(Zss[ktrace],decimate)    
                        Eds[ktrace]=stdecimate(Eds[ktrace],decimate)
                        Nds[ktrace]=stdecimate(Nds[ktrace],decimate)
                        Zds[ktrace]=stdecimate(Zds[ktrace],decimate)
                        
                    
                    ktrace+=1            
        #Read time series
        for ksta in range(Nsta):
            edata=read(datafiles[ksta]+'.e')
            ndata=read(datafiles[ksta]+'.n')
            udata=read(datafiles[ksta]+'.u')
            if decimate!=None:
                edata[0]=stdecimate(edata[0],decimate)
                ndata[0]=stdecimate(ndata[0],decimate)
                udata[0]=stdecimate(udata[0],decimate)
            if ksta==0:
                Edata=edata.copy()

es[0].data*=100
            us[0].data*=100
           
        
        if lowpass!=None:
            print('Lowpassing')
            fsample=1./e[0].stats.delta
            e[0].data=lfilter(e[0].data,lowpass,fsample,2)
            n[0].data=lfilter(n[0].data,lowpass,fsample,2)
            u[0].data=lfilter(u[0].data,lowpass,fsample,2)
            es[0].data=lfilter(es[0].data,lowpass,fsample,2)
            ns[0].data=lfilter(ns[0].data,lowpass,fsample,2)
            us[0].data=lfilter(us[0].data,lowpass,fsample,2)
        if decimate!=None:
            n[0]=stdecimate(n[0],decimate)
            e[0]=stdecimate(e[0],decimate)
            u[0]=stdecimate(u[0],decimate)

        if scale!=None:
            n[0].data=n[0].data/scale
            ns[0].data=ns[0].data/scale
            e[0].data=e[0].data/scale
            es[0].data=es[0].data/scale
            u[0].data=u[0].data/scale
            us[0].data=us[0].data/scale
        #Make plot
        if nsta>1:
            axn=axarr[k,0]
            axe=axarr[k,1]
            axu=axarr[k,2]
        else:
            axn=axarr[0]

suffix=synthsuffix
    i=where(gf[:,kgf]==1)[0]
    for k in range(len(i)):
        print('Working on '+sta[i[k]])
        if d_or_s.lower()=='d': #Read data
            n=read(path+sta[i[k]]+'.'+suffix+'.n')
            e=read(path+sta[i[k]]+'.'+suffix+'.e')
            u=read(path+sta[i[k]]+'.'+suffix+'.u')
            outname=sta[i[k]]+'.'+suffix+'.psd'
            if lowpass!=None:
                fsample=1./e[0].stats.delta
                e[0].data=lfilter(e[0].data,lowpass,fsample,10)
                n[0].data=lfilter(n[0].data,lowpass,fsample,10)
                u[0].data=lfilter(u[0].data,lowpass,fsample,10)
            if decimate!=None:
                n[0]=stdecimate(n[0],decimate)
                e[0]=stdecimate(e[0],decimate)
                u[0]=stdecimate(u[0],decimate)
        else: #Read synthetics
            n=read(path+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+suffix+'.n.sac')
            e=read(path+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+suffix+'.e.sac')
            u=read(path+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+suffix+'.u.sac')
            outname=run_name+'.'+run_number+'.'+sta[i[k]]+'.'+suffix+'.psd'
        #Compute spectra
        fn, npsd, nu = tsa.multi_taper_psd(n[0].data,Fs=1./n[0].stats.delta,adaptive=True,jackknife=False,low_bias=True)
        fe, epsd, nu = tsa.multi_taper_psd(e[0].data,Fs=1./e[0].stats.delta,adaptive=True,jackknife=False,low_bias=True)
        fu, upsd, nu = tsa.multi_taper_psd(u[0].data,Fs=1./u[0].stats.delta,adaptive=True,jackknife=False,low_bias=True)
        #Convert to dB
        npsd=10*log10(npsd)
        epsd=10*log10(epsd)
        upsd=10*log10(upsd)
        #Write to file

for ksta in range(len(i)):
        if quiet==False:  
            print('Assembling acceleration waveforms from '+stations[i[ksta]]+' into data vector.')
        n=read(GFfiles[i[ksta],kgf]+'.n')
        e=read(GFfiles[i[ksta],kgf]+'.e')
        u=read(GFfiles[i[ksta],kgf]+'.u')
        if velocity_bandpass is not None: #Apply low pass filter to data
            fsample=1./n[0].stats.delta
            n[0].data=lfilt(n[0].data,velocity_bandpass,fsample,2)
            e[0].data=lfilt(e[0].data,velocity_bandpass,fsample,2)
            u[0].data=lfilt(u[0].data,velocity_bandpass,fsample,2)
        #Decimate
        if decimate!=None:
            n[0]=stdecimate(n[0],decimate)
            e[0]=stdecimate(e[0],decimate)
            u[0]=stdecimate(u[0],decimate)
        #Make sure they are rounded to a dt interval
        dt=e[0].stats.delta
        e[0].stats.starttime=round_time(e[0].stats.starttime,dt)
        n[0].stats.starttime=round_time(n[0].stats.starttime,dt)
        u[0].stats.starttime=round_time(u[0].stats.starttime,dt)
        dvel=append(dvel,r_[n[0].data,e[0].data,u[0].data])
    
    #Tsunami
    kgf=3
    dtsun=array([])
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        if quiet==False:  
            print('Assembling tsunami waveforms from '+stations[i[ksta]]+' into data vector.')
        tsun=read(GFfiles[i[ksta],kgf])
        if tsunami_bandpass is not None: #Apply low pass filter to data

#Get RMS (L2 norm) of weigtehd data misfit
        L2static=norm(wds[kstart:kend]-wd[kstart:kend])
   
    #Displacement
    kstart=kend
    kgf=1
    i=where(GF[:,kgf]==1)[0]
    VRdisp=nan
    L2disp=nan
    if len(i)>0:
        for ksta in range(len(i)):
            sta=stations[i[ksta]]
            n=read(GFfiles[i[ksta],kgf]+'.n')
            if decimate != None:
                n[0]=stdecimate(n[0],decimate)
            npts=n[0].stats.npts
            kend+=3*npts
        #Variance reduction
        res=((d[kstart:kend]-ds[kstart:kend])**2)**0.5
        dnorm=(d[kstart:kend]**2)**0.5 #Yes i know this is dumb, shush
        VRdisp=(1-(res.sum()/dnorm.sum()))*100
        
        #Get RMS (L2 norm) of weigtehd data misfit
        L2disp=norm(wds[kstart:kend]-wd[kstart:kend])
        
    #Velocity
    kstart=kend
    kgf=2
    i=where(GF[:,kgf]==1)[0]
    VRvel=nan
    L2vel=nan

sigman=r_[sigman,ones(Nt)*sn[i[k]]]
                sigmae=r_[sigmae,ones(Nt)*se[i[k]]]
                sigmau=r_[sigmau,ones(Nt)*su[i[k]]]
        sigma_disp=r_[sigman,sigmae,sigmau]
    #Get velocioty covariance
    sigma_vel=[]
    i=where(gflist[:,2]==1)[0]
    if len(i)>0:
        #Read variances
        sn=genfromtxt(gf_file,usecols=19)
        se=genfromtxt(gf_file,usecols=20)
        su=genfromtxt(gf_file,usecols=21)
        for k in range(len(i)):
            #Get length of time series
            st=read(data_files[i[k],1]+'.n')
            st[0]=stdecimate(st[0],decimate)
            Nt=st[0].stats.npts
            #Cocnatenate
            if k==0:
                sigman=ones(Nt)*sn[i[k]]
                sigmae=ones(Nt)*se[i[k]]
                sigmau=ones(Nt)*su[i[k]]
            else:
                sigman=r_[sigman,ones(Nt)*sn[i[k]]]
                sigmae=r_[sigmae,ones(Nt)*se[i[k]]]
                sigmau=r_[sigmau,ones(Nt)*su[i[k]]]
        sigma_vel=r_[sigman,sigmae,sigmau]
    Cd=diag(r_[sigma_static,sigma_disp,sigma_vel])
    return Cd

#Static weights
    kgf=0
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        w[kinsert]=1/weights[i[ksta],0] #North
        w[kinsert+1]=1/weights[i[ksta],1] #East
        w[kinsert+2]=1/weights[i[ksta],2] #Up
        kinsert=kinsert+3
    #Displacement waveform weights
    kgf=1
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        #Read waveform to determine length of insert
        st=read(GFfiles[i[ksta],kgf]+'.n')
        if decimate!=None:
            st[0]=stdecimate(st[0],decimate)
        nsamples=st[0].stats.npts
        wn=(1/weights[i[ksta],3])*ones(nsamples)
        w[kinsert:kinsert+nsamples]=wn
        kinsert=kinsert+nsamples
        we=(1/weights[i[ksta],4])*ones(nsamples)
        w[kinsert:kinsert+nsamples]=we
        kinsert=kinsert+nsamples
        wu=(1/weights[i[ksta],5])*ones(nsamples)
        w[kinsert:kinsert+nsamples]=wu
        kinsert=kinsert+nsamples
    #velocity waveform weights
    kgf=2
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        #Read waveform to determine length of insert
        st=read(GFfiles[i[ksta],kgf]+'.n')

i=i[j]
    nsta=len(i)
    fig, axarr = plt.subplots(nsta, 3)  
    for k in range(len(i)):
        n=read(datapath+sta[i[k]]+'.'+datasuffix+'.n')
        e=read(datapath+sta[i[k]]+'.'+datasuffix+'.e')
        u=read(datapath+sta[i[k]]+'.'+datasuffix+'.u')
        if lowpass!=None:
            fsample=1./e[0].stats.delta
            e[0].data=lfilter(e[0].data,lowpass,fsample,2)
            n[0].data=lfilter(n[0].data,lowpass,fsample,2)
            u[0].data=lfilter(u[0].data,lowpass,fsample,2)
        if decimate!=None:
            n[0]=stdecimate(n[0],decimate)
            e[0]=stdecimate(e[0],decimate)
            u[0]=stdecimate(u[0],decimate)
        if scale!=None:
            n[0].data=n[0].data/scale
            e[0].data=e[0].data/scale
            u[0].data=u[0].data/scale
        #Make plot
        if nsta>1:
            axn=axarr[k,0]
            axe=axarr[k,1]
            axu=axarr[k,2]
        else:
            axn=axarr[0]
            axe=axarr[1]
            axu=axarr[2]
        axn.plot(n[0].times(),n[0].data,'k')
        axn.grid(which='both')
        axe.plot(e[0].times(),e[0].data,'k')

n.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.disp.n.sac',format='SAC')
            e.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.disp.e.sac',format='SAC')
            u.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.disp.u.sac',format='SAC')
    #Velocity
    kgf=2
    i=where(GF[:,kgf]==1)[0]
    if len(i)>0:
        for ksta in range(len(i)):
            sta=stations[i[ksta]]
            n=read(GFfiles[i[ksta],kgf]+'.n')
            e=read(GFfiles[i[ksta],kgf]+'.e')
            u=read(GFfiles[i[ksta],kgf]+'.u')
            if decimate !=None:
                n[0]=stdecimate(n[0],decimate)
                e[0]=stdecimate(e[0],decimate)
                u[0]=stdecimate(u[0],decimate)
            npts=n[0].stats.npts
            n[0].data=squeeze(ds[kinsert:kinsert+npts])
            e[0].data=squeeze(ds[kinsert+npts:kinsert+2*npts])
            u[0].data=squeeze(ds[kinsert+2*npts:kinsert+3*npts])
            kinsert+=3*npts
            n.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.vel.n.sac',format='SAC')
            e.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.vel.e.sac',format='SAC')
            u.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.vel.u.sac',format='SAC')
    #Tsunami
    kgf=3
    i=where(GF[:,kgf]==1)[0]
    if len(i)>0:
        for ksta in range(len(i)):
            sta=stations[i[ksta]]
            tsun=read(GFfiles[i[ksta],kgf])
            npts=tsun[0].stats.npts