-
Notifications
You must be signed in to change notification settings - Fork 18
/
ArrayFuncs.jl
110 lines (95 loc) · 3.65 KB
/
ArrayFuncs.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
import SeisIO: demean, demean!, taper, taper!,detrend, detrend!
export detrend, detrend!, taper, taper!, demean, demean!, hanningwindow
# Signal processing functions for arrays (rather than SeisData or SeisChannel)
"""
detrend!(X::AbstractArray{<:AbstractFloat})
Remove linear trend from array `X` using least-squares regression.
"""
function detrend!(X::AbstractArray{<:AbstractFloat})
T = eltype(X)
N = size(X,1)
# create linear trend matrix
A = similar(X,T,N,2)
A[:,2] .= T(1)
A[:,1] .= range(T(0),T(1),length=N)
# create linear trend matrix
R = transpose(A) * A
# do the matrix inverse for 2x2 matrix
# this is really slow on GPU
Rinv = inv(Array(R)) |> typeof(R)
factor = Rinv * transpose(A)
# remove trend
X .-= A * (factor * X)
return nothing
end
detrend(A::AbstractArray{<:AbstractFloat}) = (U = deepcopy(A);
detrend!(U);return U)
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)
"""
demean!(A::AbstractArray{<:AbstractFloat})
Remove mean from array `A`.
"""
function demean!(A::AbstractArray{<:AbstractFloat}; dims=1)
μ = mean(A,dims=dims)
A .-= μ
return nothing
end
demean(A::AbstractArray{<:AbstractFloat}; dims=1) = (U = deepcopy(A);
demean!(U,dims=dims);return U)
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)
"""
taper!(A,fs; max_percentage=0.05, max_length=20.)
Taper a time series `A` with sampling_rate `fs`.
Defaults to 'hann' window. Uses smallest of `max_percentage` * `fs`
or `max_length`.
# Arguments
- `A::AbstractArray`: Time series.
- `fs::Float64`: Sampling rate of time series `A`.
- `max_percentage::float`: Decimal percentage of taper at one end (ranging
from 0. to 0.5).
- `max_length::Float64`: Length of taper at one end in seconds.
"""
function taper!(A::AbstractArray{<:AbstractFloat}, fs::Float64;
max_percentage::Float64=0.05, max_length::Float64=20.)
Nrows = size(A,1)
T = eltype(A)
wlen = min(Int(floor(Nrows * max_percentage)), Int(floor(max_length * fs)), Int(
floor(Nrows/2)))
taper_sides = hanningwindow(A,2 * wlen)
A[1:wlen,:] .*= taper_sides[1:wlen]
A[end-wlen:end,:] .*= taper_sides[wlen:end]
return nothing
end
taper(A::AbstractArray{<:AbstractFloat}, fs::Float64;
max_percentage::Float64=0.05,max_length::Float64=20.) = (U = deepcopy(A);
taper!(U,fs,max_percentage=max_percentage,max_length=max_length);return U)
taper!(R::RawData; max_percentage::Float64=0.05,
max_length::Float64=20.) = taper!(R.x,R.fs,max_percentage=max_percentage,
max_length=max_length)
taper(R::RawData; max_percentage::Float64=0.05,
max_length::Float64=20.) = (U = deepcopy(R); taper!(U.x,U.fs,
max_percentage=max_percentage,max_length=max_length); return U)
taper!(C::CorrData; max_percentage::Float64=0.05,
max_length::Float64=20.) = taper!(C.corr,C.fs,max_percentage=max_percentage,
max_length=max_length)
taper(C::CorrData; max_percentage::Float64=0.05,
max_length::Float64=20.) = (U = deepcopy(C); taper!(U.corr,U.fs,
max_percentage=max_percentage,max_length=max_length); return U)
"""
hanningwindow(A,n)
Generate hanning window of length `n`.
Hanning window is sin(n*pi)^2.
"""
function hanningwindow(A::AbstractArray, n::Int)
T = eltype(A)
win = similar(A,T,n)
win .= T(pi) .* range(T(0.),stop=T(n),length=n)
win .= sin.(win).^2
return win
end