Add NodalGPU to dev (#56)

* Update Nodal stuff

* Fix CorrData filtering typo

src/ArrayFuncs.jl

@@ -33,6 +33,8 @@ detrend!(R::RawData) = detrend!(R.x)
 detrend(R::RawData) = (U = deepcopy(R); detrend!(U.x); return U)
 detrend!(C::CorrData) = detrend!(C.corr)
 detrend(C::CorrData) = (U = deepcopy(C); detrend!(U.corr); return U)
+detrend!(N::NodalData) = detrend!(N.data)
+detrend(N::NodalData) = (U = deepcopy(N); detrend!(U.data); return U)
 
 
 """
@@ -51,6 +53,8 @@ demean!(R::RawData) = demean!(R.x)
 demean(R::RawData) = (U = deepcopy(R); demean!(U.x); return U)
 demean!(C::CorrData) = demean!(C.corr)
 demean(C::CorrData) = (U = deepcopy(C); demean!(U.corr); return U)
+demean!(N::NodalData) = demean!(N.data)
+demean(N::NodalData) = (U = deepcopy(N); demean!(U.data); return U)
 
 """
    taper!(A,fs; max_percentage=0.05, max_length=20.)
@@ -92,6 +96,12 @@ taper!(C::CorrData; max_percentage::AbstractFloat=0.05,
 taper(C::CorrData; max_percentage::AbstractFloat=0.05,
       max_length::Real=20.) = (U = deepcopy(C); taper!(U.corr,U.fs,
       max_percentage=max_percentage,max_length=max_length); return U)
+taper!(N::NodalData; max_percentage::AbstractFloat=0.05,
+      max_length::Real=20.) = taper!(N.data,N.fs[1],max_percentage=max_percentage,
+      max_length=max_length)
+taper(N::NodalData; max_percentage::AbstractFloat=0.05,
+      max_length::Real=20.) = (U = deepcopy(N); taper!(U.data,U.fs[1],
+      max_percentage=max_percentage,max_length=max_length); return U)
 
 
 """

src/SeisNoise.jl

@@ -1,6 +1,6 @@
 module SeisNoise
 
-using Dates, DataFrames, DSP, FFTW, Glob, JLD2, LinearAlgebra, SeisIO
+using Dates, DataFrames, DSP, FFTW, Glob, JLD2, LinearAlgebra, SeisIO, SeisIO.Nodal
 using Statistics, StatsBase, Plots, Distributed
 using CUDA, Adapt, GPUArrays
 
@@ -13,6 +13,7 @@ end
 
 # import types first
 include("Types/NoiseData.jl")
+include("Types/NodalData.jl")
 include("Types/show.jl")
 
 # import pre and post processing tools
@@ -28,6 +29,7 @@ include("stacking.jl")
 # import  routines for doin' stuff
 include("compute_fft.jl")
 include("correlation.jl")
+include("nodalcorrelation.jl")
 include("rotation.jl")
 include("Plotting/plotting.jl")
 

src/Types/NodalData.jl

@@ -0,0 +1,80 @@
+export gpu,cpu, NodalCorr
+
+mutable struct NodalCorr <: NoiseData
+    n::Int64
+    id::Array{String,1}                 # id
+    name::Array{String,1}               # name
+    loc::Array{<:InstrumentPosition,1}  # loc
+    fs::Float64                         # fs
+    src::String                         # src
+    t::Array{Array{Int64,2},1}          # time
+    corr::AbstractArray{Float32, 2}     # actual data
+
+    function NodalCorr(
+        n        ::Int64,
+        id       ::Array{String,1},
+        name     ::Array{String,1},
+        loc      ::Array{<:InstrumentPosition,1},
+        fs       ::Float64,
+        src      ::String,
+        t        ::Array{Array{Int64,2},1} ,
+        corr     ::AbstractArray{<:AbstractFloat,2},
+        )
+
+        return new(n, id, name, loc, fs, src, t, corr)
+    end
+end
+
+NodalCorr(;
+          n        ::Int64                       = 0,
+          id       ::Array{String,1}             = Array{String}(undef,0),
+          name     ::Array{String,1}             = Array{String}(undef,0),
+          loc      ::Array{InstrumentPosition,1} = Array{InstrumentPosition}(undef,0), 
+          fs       ::Float64                     = zero(Float64),
+          src      ::String                      = "",
+          t        ::Array{Array{Int64,2},1}     = Array{Array{Int64,2},1}(undef,0),
+          corr     ::AbstractArray{<:AbstractFloat,2} = Array{Float32,2}(undef, 0, 2),
+          ) = CorrData(n, id, name, loc, fs, src, t, corr)
+
+
+# functions for adapting NodalData to GPU 
+cpu(m::NodalData) = adapt(Array,m)
+gpu(x::NodalData) = use_cuda[] ? cu(x) : x
+
+# function for adapting NodalCorr to GPU 
+cpu(m::NodalCorr) = adapt(Array,m)
+gpu(x::NodalCorr) = use_cuda[] ? cu(x) : x
+
+Adapt.adapt_structure(to, N::NodalData) = NodalData(
+    N.n,
+    N.ox,
+    N.oy,
+    N.oz,
+    N.info,
+    N.id,
+    N.name,
+    N.loc,
+    N.fs,
+    N.gain,
+    N.resp,
+    N.units,
+    N.src,
+    N.misc,
+    N.notes,
+    N.t,
+    adapt(to,N.data),
+)
+
+Adapt.adapt_structure(to, N::NodalCorr) = NodalCorr(
+    N.n,
+    N.id,
+    N.name,
+    N.loc,
+    N.fs,
+    N.gain,
+    N.resp,
+    N.units,
+    N.src,
+    N.t,
+    adapt(to,N.corr),
+)

src/correlation.jl

@@ -268,6 +268,18 @@ end
 whiten(R::RawData,freqmin::Real, freqmax::Real; pad::Int=50) =
       (U = deepcopy(R); whiten!(U,freqmin,freqmax,pad=pad);return U)
 
+function whiten!(N::NodalData,freqmin::Real, freqmax::Real; pad::Int=50)
+    @assert freqmin > 0 "Whitening frequency must be greater than zero."
+    @assert freqmax <= N.fs[1] / 2 "Whitening frequency must be less than or equal to Nyquist frequency."
+    Npts = size(N.data,1) # number of data points
+    FFT = rfft(N.data,1)
+    whiten!(FFT,freqmin,freqmax,N.fs[1], Npts, pad=pad)
+    N.data .= irfft(FFT,Npts,1)
+    return nothing
+end
+whiten(N::NodalData,freqmin::Real, freqmax::Real; pad::Int=50) =
+    (U = deepcopy(N); whiten!(U,freqmin,freqmax,pad=pad);return U)
+
 """
 
   coherence!(F,half_win, water_level)

src/distance.jl

@@ -374,3 +374,22 @@ function get_loc(C::CorrData)
 
     return C.loc, loc2
 end
+
+function get_azi(loc1::NodalLoc,loc2::NodalLoc)
+    azi = atan(loc2.y - loc1.y, loc2.x - loc1.x) * 180 / π 
+    if azi <= 90 
+        azi = 90 - azi 
+    else
+        azi = 450 - azi 
+    end
+    return azi  
+end
+
+function get_azi(N::NodalData)
+    out = zeros(N.n)
+    for ii = 1:N.n-1
+        out[ii] = get_azi(N.loc[ii],N.loc[ii+1])
+    end
+    out[end] = out[end-1]
+    return out
+end

src/filter.jl

@@ -72,9 +72,9 @@ bandpass(R::RawData,freqmin::Real,freqmax::Real;
         setfield!(U,:freqmin,Float64(freqmin));
         setfield!(U,:freqmax,Float64(min(freqmax,U.fs/2)));return U)
 bandpass!(C::CorrData,freqmin::Real,freqmax::Real;
-          corners::Int=4,zerophase::Bool=true) = (bandpass!(C.corr,freqmin,freqmax,
-          C.fs,corners=corners,zerophase=zerophase);setfield!(C,:freqmin,Float64(freqmin));
-          setfield!(C,:freqmax,Float64(min(freqmax,C.fs/2)));return nothing)
+        corners::Int=4,zerophase::Bool=true) = (bandpass!(C.corr,freqmin,freqmax,
+        C.fs,corners=corners,zerophase=zerophase);setfield!(C,:freqmin,Float64(freqmin));
+        setfield!(C,:freqmax,Float64(min(freqmax,C.fs/2)));return nothing)
 bandpass(C::CorrData,freqmin::Real,freqmax::Real;
         corners::Int=4,zerophase::Bool=true) = (U = deepcopy(C);bandpass!(U.corr,
         freqmin,freqmax,U.fs,corners=corners,zerophase=zerophase);
@@ -91,6 +91,13 @@ bandpass!(S::SeisData, freqmin::Real, freqmax::Real;
                    fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandpass")
 bandpass(S::SeisData, freqmin::Real, freqmax::Real; corners::Int=4,
          zerophase::Bool=true) = filtfilt(S,fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandpass")
+bandpass!(N::NodalData,freqmin::Real,freqmax::Real;
+         corners::Int=4,zerophase::Bool=true) = (bandpass!(N.data,freqmin,freqmax,
+         N.fs[1],corners=corners,zerophase=zerophase);return nothing)
+bandpass(N::NodalData,freqmin::Real,freqmax::Real;
+       corners::Int=4,zerophase::Bool=true) = (U = deepcopy(N);bandpass!(N.data,
+       freqmin,freqmax,N.fs[1],corners=corners,zerophase=zerophase);
+       return U)
 
 """
   bandstop!(A,freqmin,freqmax,fs,corners=4,zerophase=true)
@@ -175,6 +182,12 @@ bandstop!(S::SeisData, freqmin::Real, freqmax::Real;
 bandstop(S::SeisData,freqmin::Real, freqmax::Real; corners::Int=4,
          zerophase::Bool=true) = filtfilt(S,
          fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandstop")
+bandstop!(N::NodalData,freqmin::Real,freqmax::Real;
+    corners::Int=4,zerophase::Bool=true) = bandstop!(N.data,freqmin,freqmax,
+    N.fs[1],corners=corners,zerophase=zerophase)
+bandstop(N::NodalData,freqmin::Real,freqmax::Real;
+    corners::Int=4,zerophase::Bool=true) = (U = deepcopy(N);bandstop!(N.data,
+    freqmin,freqmax,U.fs[1],corners=corners,zerophase=zerophase);return U)
 
 """
 lowpass(A,freq,fs,corners=4,zerophase=true)
@@ -250,6 +263,12 @@ lowpass!(S::SeisData, freq::Real;corners::Int=4, zerophase::Bool=true) =
     filtfilt!(S,fh=Float64(freq),np=corners,rt="Lowpass")
 lowpass(S::SeisData,freq::Real; corners::Int=4,zerophase::Bool=true) =
     filtfilt(S,fh=Float64(freq),np=corners,rt="Lowpass")
+lowpass!(N::NodalData,freq::Real; corners::Int=4,
+    zerophase::Bool=true) = (lowpass!(N.data,freq,N.fs[1],corners=corners,
+    zerophase=zerophase);return nothing)
+lowpass(N::NodalData,freq::Real; corners::Int=4,
+    zerophase::Bool=true) = (U = deepcopy(N);lowpass!(N.data,freq,U.fs[1],
+    corners=corners,zerophase=zerophase);return U)
 
 """
    highpass(A,freq,fs,corners=4,zerophase=true)
@@ -321,6 +340,12 @@ highpass!(S::SeisData, freq::Real;corners::Int=4, zerophase::Bool=true) =
     filtfilt!(S,fl=Float64(freq),np=corners,rt="Highpass")
 highpass(S::SeisData,freq::Real; corners::Int=4,zerophase::Bool=true) =
     filtfilt(S,fl=Float64(freq),np=corners,rt="Highpass")
+highpass!(N::NodalData,freq::Real; corners::Int=4,
+    zerophase::Bool=true) = (highpass!(N.data,freq,N.fs[1],corners=corners,
+    zerophase=zerophase);return nothing)
+highpass(N::NodalData,freq::Real; corners::Int=4,
+    zerophase::Bool=true) = (U = deepcopy(N);highpass!(N.data,freq,U.fs[1],
+    corners=corners,zerophase=zerophase);return U)
 
 """
     envelope(A)

src/nodalcorrelation.jl

@@ -0,0 +1,79 @@
+export scalablexcorr, nodalxcorr, autocorrelate
+
+function scalablexcorr(A::AbstractArray,Ntau::Int,kthresh::AbstractFloat)
+    V,S,U = svd(A)
+    kind = findall(S .> maximum(S) * 0.05)
+    kmax = maximum(kind)
+
+    Nt,Ns = size(A)
+
+    # check that Ntau < Nt / 2 
+    Ntau = min(Ntau,Nt ÷ 2)
+
+    # preallocate view 
+    V0 = @view V[Ntau÷2:Nt-Ntau÷2-1,1:kmax]
+    U0 = @view U[:,1:kmax]
+    if isa(V,AbstractGPUArray)
+        VT0 = transpose(V[Ntau÷2:Nt-Ntau÷2-1,1:kmax])
+        UT0 = transpose(U[:,1:kmax])
+    else
+        VT0 = V0'
+        UT0 = U0'
+    end
+
+    Cout = similar(A,Ns,Ns,Ntau)
+
+    for itau = -Ntau ÷ 2 + 1 : Ntau ÷ 2 
+        W = VT0 * V[itau + Ntau÷2:Nt + itau - Ntau÷2-1,1:kmax]
+        Cout[:,:,itau+Ntau÷2] .= U0 * W * UT0
+    end
+    return Cout
+end
+
+function scalablexcorr(N::NodalData,Ntau::Real;kthresh=0.05)
+    if isa(Ntau,AbstractFloat)
+        Ntau = convert(Int,round(Ntau * N.fs[1]))
+    end
+    return scalablexcorr(N.data,Ntau,kthresh)
+end
+
+function nodalxcorr(N::NodalData,maxlag::Real)
+    if isa(maxlag,AbstractFloat)
+        maxlag = convert(Int,round(maxlag * N.fs[1]))
+    end
+    return nodalxcorr(N.data,maxlag)
+
+end
+
+function nodalxcorr(A::AbstractArray,maxlag::Int)
+    Nt,Ns = size(A)
+    Ncorr = Ns * (Ns -1) ÷ 2 
+    Cout = similar(A,maxlag * 2 + 1,Ncorr)
+    FFTs = rfft(A,1)
+    cstart = 0
+    cend = 0
+    for ii = 1:Ns-1
+            cstart = cend + 1 
+            cend = cstart + Ns - ii - 1
+            Cout[:,cstart:cend] .= correlate(FFTs[:,ii],FFTs[:,ii+1:end],Nt,maxlag)
+    end
+    return Cout
+end
+
+function autocorrelate(N::NodalData,maxlag::Real)
+    if isa(maxlag,AbstractFloat)
+        maxlag = convert(Int,round(maxlag * N.fs[1]))
+    end
+    return autocorrelate(N.data,maxlag)
+end
+
+function autocorrelate(A::AbstractArray,maxlag::Int64)
+    Nt,Ns = size(A)
+    FFT = rfft(A,1)
+    corrT = irfft(conj.(FFT) .* FFT,Nt,1)
+    # return corr[-maxlag:maxlag]
+    t = vcat(0:Int(Nt  / 2)-1, -Int(Nt  / 2):-1)
+    ind = findall(abs.(t) .<= maxlag)
+    newind = fftshift(ind,1)
+    return corrT[newind,:]
+end