-
Notifications
You must be signed in to change notification settings - Fork 18
/
phase_shift.jl
76 lines (62 loc) · 2.28 KB
/
phase_shift.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
export phase_shift, phase_shift!
"""
phase_shift!(C::SeisChannel)
Phase shift SeisChannel if starttime is not aligned with sampling period.
For example, a channel with starttime of 2016-07-04T00:00:00.008 and a sampling
rate of 100 Hz, the starttime is off the sampling period (0.01) by 2 millisecond
and should be at 2016-07-04T00:00:00.01. The waveform will be phase shifted by
2 milliseconds in the frequency domain. If ϕshift = false, will only correct
starttime of channel.
# Arguments
- `C::SeisChannel`: SeisChannel.
- `ϕshift::Bool`: Boolean to allow/disable phase-shifting.
"""
function phase_shift!(C::SeisChannel; ϕshift::Bool=true)
# get time offset from sampling rate
dt = Float32(1. / C.fs)
off = mod(mod(C.t[1,2],1000000)*Float32(1e-6),dt)
n = length(C.x)
if dt - off <= eps(eltype(off))
off = 0.
end
if off != 0.
if off <= dt / 2.
off = -off
else
off = dt - off
end
# phase shift data
if ϕshift
freq = FFTW.rfftfreq(length(C.x),C.fs)
fftdata = rfft(C.x)
fftdata .= fftdata .* exp.(1im .* 2 .* pi .* freq .* off)
C.x[:] = irfft(fftdata,length(C.x))
end
# reset time of signal
C.t[1,2] += off * 1e6
end
return nothing
end
phase_shift(C::SeisChannel;ϕshift::Bool=true) = (U = deepcopy(C);
phase_shift!(U,ϕshift=ϕshift);return U)
"""
phase_shift!(S::SeisData)
For example, a channel with starttime of 2016-07-04T00:00:00.008 and a sampling
rate of 100 Hz, the starttime is off the sampling period (0.01) by 2 millisecond
and should be at 2016-07-04T00:00:00.01. The waveform will be phase shifted by
2 milliseconds in the frequency domain. If ϕshift = false, will only correct
starttime of channel.
This function will not work properly for sampling rates >= 1000 Hz due to Julia
DateTime resolving time only to milliseconds!
# Arguments
- `S::SeisData`: SeisData.
- `ϕshift::Bool`: Boolean to allow/disable phase-shifting.
"""
function phase_shift!(S::SeisData;ϕshift::Bool=true)
@inbounds for i = 1:S.n
phase_shift!(S[i],ϕshift=ϕshift)
end
return nothing
end
phase_shift(S::SeisData;ϕshift::Bool=true) = (U = deepcopy(S);
phase_shift!(U,ϕshift=ϕshift);return U)