-
Notifications
You must be signed in to change notification settings - Fork 18
/
tools.jl
106 lines (91 loc) · 2.96 KB
/
tools.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
export snr, smooth, smooth!, abs_max, abs_max!, standardize, standardize!
export mad, std_threshold
"""
snr(A)
Signal to noise ratio of cross-correlations in matrix `A`.
Follows method of Clarke et. al, 2011. Measures SNR at each point.
"""
function snr(A::AbstractArray)
Nrows,Ncols = size(A)
A_mean = mean(A,dims=2)
# calculate noise and envelope functions
sigma = mean(A.^2,dims=2) .- A_mean.^2
sigma .= sqrt.(sigma./ (Ncols-1))
s = abs.(A_mean .+ im .* hilbert(A_mean))
return s ./ sigma
end
snr(C::CorrData) = snr(C.corr)
"""
smooth(A, half_win)
Smooth array `A` with half-window `half_win` (defaults to 3).
"""
function smooth!(A::AbstractArray, half_win::Int=3, dims::Int=1)
T = eltype(A)
window_len = 2 * half_win + 1
csumsize = tuple(collect(size(A)) .+ [i==1 ? 2 * (window_len - 1) + 1 : 0 for i in 1:ndims(A)]...)
csum = similar(A,T,csumsize)
csum[1:window_len,:] .= zero(T)
csum[end - window_len + 1:end,:] .= zero(T)
csum[window_len+1:end-window_len + 1,:] .= A
csum .= cumsum(csum,dims=dims)
weight = similar(A,T,size(A,1))
weight[1:half_win] = T.(window_len ÷ 2 + 1:window_len - 1)
weight[half_win + 1: end - half_win] .= T(window_len)
weight[end-half_win:end] = T.(window_len:-1:window_len ÷ 2 + 1)
A[:,:] .= (csum[window_len+half_win+1:end-half_win,:] .- csum[half_win+1:end-window_len-half_win,:]) ./ weight
return nothing
end
smooth(A::AbstractArray,half_win::Int=3, dims::Int=1) =
(U = deepcopy(A);smooth!(U,half_win,dims);return U)
"""
abs_max!(A)
Returns array `A` divided by its absolute maximum value.
"""
function abs_max!(A::AbstractArray)
maxabs = maximum(abs.(A),dims=1)
if any(maxabs .== 0)
throw(DomainError("All zero column leads to NaN"))
end
A ./= maxabs
return nothing
end
abs_max(A::AbstractArray) = (U = deepcopy(A);abs_max!(U);return U)
abs_max!(C::CorrData) = abs_max!(C.corr)
abs_max(C::CorrData) = (U = deepcopy(C);abs_max!(U);return U)
"""
standardize!(A)
Demean and standardize array `A` to unit std.
"""
function standardize!(A::AbstractArray)
A .-= mean(A,dims=1)
A ./= std(A,dims=1)
return nothing
end
standardize(A::AbstractArray) = (U = deepcopy(A);standardize!(U);return U)
standardize!(C::CorrData) = standardize!(C.corr)
standardize(C::CorrData) = (U = deepcopy(C);standardize!(U);return U)
"""
mad(A)
Median Absolute Deviation of array `A`.
MAD = median(|Xi- median(A)|)
"""
function mad(A::AbstractArray)
return median(abs.(A .- median(A,dims=1)),dims=1)
end
mad(C::CorrData) = mad(C.corr)
"""
std_threshold(A,thresh)
Returns indices of cols of `A` where max(abs.A[:,col]))/std(A[:,col]) < `max_std`.
"""
function std_threshold(A::AbstractArray, max_std::Real)
stds = std(A,dims=1)[1,:]
maxs = maximum(abs.(A),dims=1)[1,:]
threshs = maxs ./ stds
ind = []
for ii = 1:length(threshs)
if threshs[ii] < max_std
append!(ind,ii)
end
end
return ind
end