jpjones76 / SeisIO.jl

Julia language support for geophysical time series data
http://seisio.readthedocs.org
Other
47 stars 21 forks source link

Allow for upsampliing in resample! #50

Closed tclements closed 4 years ago

tclements commented 4 years ago

resample!(S,fs=fs throws a BoundsError when trying to upsample data where fs > S.fs.

MWE:

using SeisIO 
S = SeisIO.RandSeis.randSeisData(3,nx=2^14)
SeisData with 3 channels (3 shown)
    ID: W7.IJVJP.M4.CH3                    0M.OVX..STZ                        3A.GCPM..MLE                       
  NAME: 3wVmTMZQ2g0Z6Ls4Axh2JlB4nuSbO7lU5… G103gKlj6wY2d1sQGYSprCrwLePPGzZWC… 9LiHmHKS8LR8pe7pTHOE8pA0ZHvNFg     
   LOC: [-30.1085, 50.2216, 964.76, 67.77… [-15.6587, -85.0015, 480.462, 53.… -29.9009 N, 131.203 E, 921.742 m   
    FS: 250.0                              62.5                               2.0                                
  GAIN: 4.01747e8                          7.19092e9                          9.03546e7                          
  RESP: a0 1.0, f0 1.0, 8z, 8p             a0 1.0, f0 1.0, 8z, 8p             a0 1.0, f0 1.0, 8z, 8p             
 UNITS: m/s                                m                                  m                                  
   SRC: randSeisChannel(c=false, nx=16384… randSeisChannel(c=false, nx=16384… randSeisChannel(c=false, nx=16384… 
  MISC: 5 entries                          5 entries                          15 entries                         
 NOTES: 1 entries                          1 entries                          1 entries                          
     T: 2020-06-16T14:04:30 (5 gaps)       2020-06-16T14:04:30 (8 gaps)       2020-06-16T14:04:31 (7 gaps)       
     X: +1.194e+00                         +9.069e-01                         +2.358e+00                         
        +7.154e-01                         +1.946e+00                         -4.762e-01                         
            ...                                ...                                ...                            
        +3.401e-01                         -1.326e-01                         -3.853e-01                         
        (nx = 16384)                       (nx = 16384)                       (nx = 16384)  

resample(S,fs=100.)
ERROR: BoundsError: attempt to access 32770-element reinterpret(Float32, ::Array{Float64,1}) at index [1:819175]

DSP.jl's version of resample has no problem upsampling data where fs > S.fs.

resample(S[2].x,100. / S.fs[2])
26214-element Array{Float64,1}:
  0.9069498214529474
  1.716448479926333
  1.9163589276559438
  1.2573744598674774
  0.43120842346064603
  0.3527387116546516
  0.8284999681254623
  0.5564058433573061
 -0.8164715674304706
 -1.713997408026421
  ⋮
 -0.6770823868166385
  0.8964191087353106
  0.04654695052155877
 -1.4902814411925862
 -0.2151765050416976
  1.0116242219244012
 -0.34311131524348026
 -1.0160792803040768
  0.09180860342822253

I see that SeisIO's version of downsampling is much faster than DSP version:

using BenchmarkTools
# seisio resample
@benchmark resample(S[1],100.)
BenchmarkTools.Trial: 
  memory estimate:  527.19 KiB
  allocs estimate:  200
  --------------
  minimum time:     290.964 μs (0.00% GC)
  median time:      299.393 μs (0.00% GC)
  mean time:        327.932 μs (3.95% GC)
  maximum time:     2.741 ms (80.16% GC)
  --------------
  samples:          10000
  evals/sample:     1

# DSP resample
@benchmark resample(S[1].x,100. / S.fs[1])
BenchmarkTools.Trial: 
  memory estimate:  343.80 KiB
  allocs estimate:  62
  --------------
  minimum time:     1.484 ms (0.00% GC)
  median time:      1.556 ms (0.00% GC)
  mean time:        1.577 ms (0.50% GC)
  maximum time:     3.738 ms (50.90% GC)
  --------------
  samples:          3169
  evals/sample:     1

Does it make sense to call DSP.resample when S.fs > fs? Happy to submit a PR on this which calls out to DSP.resample when upscaling is needed.

jpjones76 commented 4 years ago

Hmm.... you could call DSP.resample, sure, but this appears to be an oversight on my part. I definitely intended to allow upsampling. I'm guessing that I need an explicit array resize check when S.fs > fs.

tclements commented 4 years ago

Ok, I wasn't sure if resample was intended to do upsampling too. DSP.resample not the best option then.

jpjones76 commented 4 years ago

Fixed on dev, will merge into master if tests pass.