Bold should be normalized such that ||b||=1 for the weight computation, as per Pavlis and Vernon (2010).
Here in L80 & L83 (and again when this is iterated) Bold becomes Bold ./ norm(Bold,2):
https://github.com/tclements/SeisNoise.jl/blob/86cdb8e7cfa47dbd8cbef9833d5ca1a613ef595f/src/stacking.jl#L80-L86
Removing L86 (and the normalization of w in later iterations) then gives the weights a real meaning, namely w.*r = cor(Bold, A[ii]). This explains the similarity (and gives the proportionality) in Chengxin's slide 32.

Also, P&V2010 compute the L1 norm in the numerator for convergence checks:
https://github.com/tclements/SeisNoise.jl/blob/86cdb8e7cfa47dbd8cbef9833d5ca1a613ef595f/src/stacking.jl#L91
The choice of convergence criteria is nonunique, so L2 in the numerator probably works just as well, albeit with modified eps.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Hi Tim!
It's the first time that I use Julia, and I get an error after installing "SeisNoise" successfully on my MacBook as follow:

It doesn't seem to have found the 'rfft' function. Do you konw the reason?

Thanks in advance

Fu

I updated after issue #11 was fixed. Now the demo in the SeisNoise front page is not working. (It was also not working when I tried adding SeisNoise#master as suggested to get the fix prior to the Julia registry update. Below is the error. SeisNoise compiles now, which it should after the issue #11 fix, but there is a new problem. Now there is a bounds error in compute_fft(). Any ideas? This is line 8 in the demo on the front page of this repository.

julia> FFT1 = compute_fft(S1,freqmin, freqmax, fs, cc_step, cc_len,
time_norm=false,to_whiten=false)
ERROR: BoundsError: attempt to access (Float32[8.9519034e-8 8.951522e-7 … -9.883005e-8 -1.2194332e-6; 1.3060095e-7 7.1878003e-7 … -1.3441387e-7 -9.5242154e-7; … ; -1.1920544e-7 6.684296e-7 … -5.1719906e-7 -4.7945866e-12; -1.9726252e-7 6.3828196e-7 … -5.381539e-7 4.4324054e-12], [1.13875245e9, 1.1387529e9, 1.13875335e9, 1.1387538e9, 1.13875425e9, 1.1387547e9, 1.13875515e9, 1.1387556e9, 1.13875605e9, 1.1387565e9 … 1.13883255e9, 1.138833e9, 1.13883345e9, 1.1388339e9, 1.13883435e9, 1.1388348e9, 1.13883525e9, 1.1388357e9, 1.13883615e9, 1.1388366e9])
at index [3]
Stacktrace:
[1] indexed_iterate at ./tuple.jl:63 [inlined]
[2] #compute_fft#141(::Bool, ::Bool, ::Float64, ::Bool, ::typeof(compute_fft), ::SeisData, ::Float64, ::Float64, ::Float64, ::Int64, ::Int64) at /Users/dmikesell/.julia/packages/SeisNoise/C5IgW/src/compute_fft.jl:260
[3] (::SeisNoise.var"#kw##compute_fft")(::NamedTuple{(:time_norm, :to_whiten),Tuple{Bool,Bool}}, ::typeof(compute_fft), ::SeisData, ::Float64, ::Float64, ::Float64, ::Int64, ::Int64) at ./none:0
[4] top-level scope at none:0

In whiten!() that the freqmin and freqmax have to be float values. When I used 1 instead of 1.0 for freqmax, it failed. Maybe we can add some flexibility in the input but convert it forcefully to float.

due to pinned dependencies in SeisIO. See jpjones76/SeisIO.jl#93 and proposed fix jpjones76/SeisIO.jl#94. Depending on how the SeisIO version is bumped after that fix, compat changes may be needed here too.

In a local instance of SeisNoise.jl that depends on SeisIO.jl with the proposed fix implemented, compilation of this package is successful.

Process_raw does not seem to handle SeisData objects with multiple channels, only SeisData objects with single channels (see screenshot). Might be a nice feature instead of converting to an array of single channel SeisData objects. Broadcasting would also be a nice feature, but more involved.

When using shorten! to subset a CorrData, if the newlag value is not a multiple of the sampling rate, the new maxlag should be set to the nearest common multiple of the sampling rate that is less than newlag e.g. C.maxlag = newlag - (newlag % (1/C.fs)).

MWE:

julia> C = load_corr(corrfile,comp);

julia> C.fs
40.0

julia> C.maxlag
20.0

julia> newlag = 16.01
16.01

julia> newlag % (1/C.fs) # newlag is off sampling period by 0.01 seconds
0.010000000000000675

julia> shorten!(C,newlag)

julia> C.maxlag
16.01

julia> C.maxlag == newlag 
true

In this example, the true value for C.maxlag should be 16.0.

Hi! ,

I am trying to slice a single channel into 70 secs windows but RawData() is showing InexactError. Is this because of some constraint on cc_len and cc_step values?

I recently updated to Julia 1.3.1 after seeing that SeisIO issues have been fixed. I tried to run the SeisNoise example on the github front page of the repository. I am now getting the following error.

WARNING: both FFTW and DSP export "rfftfreq"; uses of it in module SeisNoise must be qualified

Looks like there is a problem in these two packages and we need to modify SeisNoise to choose which we are using: FFTW or DSP. Have you seen this error yet and is there a plan to fix it?

The problem occurs on line 29 of phase_shift.jl

freq = rfftfreq(length(C.x),C.fs)

Dear Tim

I loaded 72 traces in sac format in S. The next step is to build a structure using R = RawData (S,cc_len, cc_step). I only see one trace in this step with the number of windows given by cc_len and cc_step. The question is whether the RawData function can operate on all channels found in S to convert all SeisData to RawData. Also, suppose we use the correlate function. In that case, ¿it is necessary to perform the individual correlation between pairs of stations (through a loop), or can you achieve the correlation between all pairs of stations and store it? Finally, ¿is possible to export Crosscorrelation results in SAC format?

Thanks in advance

Martín

On Julia 1.8.2 I run into dependency issues wrangling around hdf5. Plots is not tied in Project.toml, I am not used to Julia and this is a wild guess.

WARNING: could not import HDF5.Group into _hdf5_implementation
WARNING: could not import HDF5.Dataset into _hdf5_implementation
┌ Warning: Error requiring `HDF5` from `Plots`
│   exception =
│    LoadError: UndefVarError: Group not defined
│    Stacktrace:
│      [1] top-level scope
│        @ ~/.julia/packages/Plots/wutJB/src/backends/hdf5.jl:138
│      [2] include(mod::Module, _path::RelocatableFolders.Path)
│        @ Base ./Base.jl:419
│      [3] include(x::RelocatableFolders.Path)
│        @ Plots ~/.julia/packages/Plots/wutJB/src/Plots.jl:1
│      [4] top-level scope
│        @ ~/.julia/packages/Plots/wutJB/src/init.jl:68
│      [5] eval
│        @ ./boot.jl:368 [inlined]
│      [6] eval
│        @ ~/.julia/packages/Plots/wutJB/src/Plots.jl:1 [inlined]
│      [7] (::Plots.var"#360#405")()
│        @ Plots ~/.julia/packages/Requires/Z8rfN/src/require.jl:101
│      [8] macro expansion
│        @ ./timing.jl:382 [inlined]
│      [9] err(f::Any, listener::Module, modname::String, file::String, line::Any)
│        @ Requires ~/.julia/packages/Requires/Z8rfN/src/require.jl:47
│     [10] (::Plots.var"#359#404")()
│        @ Plots ~/.julia/packages/Requires/Z8rfN/src/require.jl:100
│     [11] withpath(f::Any, path::String)
│        @ Requires ~/.julia/packages/Requires/Z8rfN/src/require.jl:37
│     [12] (::Plots.var"#358#403")()
│        @ Plots ~/.julia/packages/Requires/Z8rfN/src/require.jl:99
│     [13] listenpkg(f::Any, pkg::Base.PkgId)
│        @ Requires ~/.julia/packages/Requires/Z8rfN/src/require.jl:20
│     [14] macro expansion
│        @ ~/.julia/packages/Requires/Z8rfN/src/require.jl:98 [inlined]
│     [15] __init__()
│        @ Plots ~/.julia/packages/Plots/wutJB/src/init.jl:67
│     [16] _include_from_serialized(pkg::Base.PkgId, path::String, depmods::Vector{Any})
│        @ Base ./loading.jl:831
│     [17] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String, build_id::UInt64)
│        @ Base ./loading.jl:1039
│     [18] _require(pkg::Base.PkgId)
│        @ Base ./loading.jl:1315
│     [19] _require_prelocked(uuidkey::Base.PkgId)
│        @ Base ./loading.jl:1200
│     [20] macro expansion
│        @ ./loading.jl:1180 [inlined]
│     [21] macro expansion
│        @ ./lock.jl:223 [inlined]
│     [22] require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:1144
│     [23] include(mod::Module, _path::String)
│        @ Base ./Base.jl:419
│     [24] exec_options(opts::Base.JLOptions)
│        @ Base ./client.jl:303
│     [25] _start()
│        @ Base ./client.jl:522
│    in expression starting at /home/marius/.julia/packages/Plots/wutJB/src/backends/hdf5.jl:35
└ @ Requires ~/.julia/packages/Requires/Z8rfN/src/require.jl:51

HDF5 update is held back by SeisNoise, thus causing incomp w newer versions of Plots

Error found by @xtyangpsp: whiten! accepts only ComplexF32 data and not ComplexF64. MWE:

using SeisNoise 
A = rand(ComplexF64,100,10);
whiten!(A,1.,2.,10.,50)
ERROR: MethodError: no method matching whiten!(::Array{Complex{Float64},2}, ::Float64, ::Float64, ::Float64, ::Int64)
Closest candidates are:
  whiten!(::AbstractArray{Complex{Float32},N} where N, ::Real, ::Real, ::Real, ::Int64; pad) at /home/timclements/.julia/packages/SeisNoise/Bosrn/src/correlation.jl:168
  whiten!(::FFTData, ::Real, ::Real; pad) at /home/timclements/.julia/packages/SeisNoise/Bosrn/src/correlation.jl:221
  whiten!(::RawData, ::Real, ::Real; pad) at /home/timclements/.julia/packages/SeisNoise/Bosrn/src/correlation.jl:245
  ...
Stacktrace:
 [1] top-level scope at REPL[4]:1

Need to rewrite whiten! kernel to accept Complex{T} where T <: AbstractFloat.

Hi Tim,

thanks for your work and this library! When I was working on https://github.com/pyrocko/lightguide, I wanted to run many (>1000) rffts in near real-time. A critical key to performance was the reuse of FFTWs plans (https://www.fftw.org/fftw3_doc/Using-Plans.html). As far as I understand FFTW.jl supports planing (https://juliamath.github.io/FFTW.jl/latest/fft.html#FFTW.plan_r2r).

That increased performance like 10x for me. The matrix size from RawData does not change, maybe you can leverage the same strategy here?

Best wishes,

Marius

Dear professor,
I want to know how to use the package of SeisNoise.jl and SeisIO to batch process cross-correlation function between lots of station Seisdata. I hvae got one hour station data cross-correlation function between two station.My codes are as follow.

(base) [Fangfan@localhost ~]$ cd /home/Fangfan/huanan2022/test
(base) [Fangfan@localhost test]$ julia
julia> using Pkg
julia> Pkg.build()
julia> using SeisNoise,SeisIO,Plots
julia> fs = 40. # sampling frequency in Hz
julia> freqmin,freqmax = 0.1,0.2 # min and max frequencies
julia> cc_step, cc_len = 450., 1800. # corrleation step and length in S
julia> maxlag = 60. # maximum lag time in correlation
julia> S1=read_data("sac","2022.07.28.10.14.50.YD202.BHE.sac")
julia> S2=read_data("sac","2022.07.28.10.12.12.YD204.BHE.sac")
julia> merge!(S1)
julia> merge!(S2)
julia> ungap!(S1)
julia> ungap!(S2)
julia> process_raw!(S1,fs)
julia> process_raw!(S2,fs)
julia> R = RawData.([S1,S2],cc_len,cc_step)
julia> detrend!.(R)
julia> taper!.(R)
julia> bandpass!.(R,freqmin,freqmax,zerophase=true)
julia> FFT = rfft.(R)
julia> whiten!.(FFT,freqmin,freqmax)
julia> C = correlate(FFT[1],FFT[2],maxlag)
julia> clean_up!(C,freqmin,freqmax)
julia> abs_max!(C)
julia> plot(C)
julia> savefig("2.png")

The attachment is the result of the program running

Have you had issues with the sampling frequency and process_raw!()? We noticed a new error recently and it took us a little bit to figure out why, because our code was previously working. I paste an example below which uses the example on the SeisNoise.jl page. I download some data and process using the actual sample frequency. That works fine. Then I download the data again and use a quarter of the that sampling frequency during the processing. Then I get a resize error. Any ideas? Have you seen this before? fs is supposed to Downsamples data to sampling rate fs`` according to compute_fft.jl. What are we screwing up here?

julia> fs = 40. # sampling frequency in Hz
40.0

julia> S1 = get_data("FDSN","CI.SDD..BHZ",src="SCEDC",s=s,t=t)

SeisData with 1 channels (1 shown)
ID: CI.SDD..BHZ
NAME: Saddleback
LOC: 33.5526 N, -117.662 E, 120.0 m
FS: 40.0
GAIN: 6.2652e8
RESP: a0 5.73035e8, f0 0.03, 2z, 5p
UNITS: m/s
SRC: http://service.scedc.caltech.edu/…
MISC: 4 entries
NOTES: 4 entries
T: 2019-02-03T00:00:00 (0 gaps)
X: -1.992e+03
-1.847e+03
...
-1.106e+03
(nx = 3456000)
C: 0 open, 0 total

julia> process_raw!(S1,fs)

julia> fs = 10.0
10.0

julia> S1 = get_data("FDSN","CI.SDD..BHZ",src="SCEDC",s=s,t=t)
SeisData with 1 channels (1 shown)
ID: CI.SDD..BHZ
NAME: Saddleback
LOC: 33.5526 N, -117.662 E, 120.0 m
FS: 40.0
GAIN: 6.2652e8
RESP: a0 5.73035e8, f0 0.03, 2z, 5p
UNITS: m/s
SRC: http://service.scedc.caltech.edu/…
MISC: 4 entries
NOTES: 4 entries
T: 2019-02-03T00:00:00 (0 gaps)
X: -1.992e+03
-1.847e+03
...
-1.106e+03
(nx = 3456000)
C: 0 open, 0 total

julia> process_raw!(S1,fs)
ERROR: cannot resize array with shared data
Stacktrace:
[1] deleteat!(::Array{Float32,1}, ::UnitRange{Int64}) at ./array.jl:838
[2] #resample!#157(::Int64, ::Float64, ::typeof(resample!), ::SeisData) at /Users/dmikesell/.julia/packages/SeisIO/mlqvk/src/Processing/resample.jl:125
[3] #resample! at ./none:0 [inlined]
[4] #process_raw!#127(::Bool, ::typeof(process_raw!), ::SeisData, ::Float64) at /Users/dmikesell/.julia/packages/SeisNoise/I9Knu/src/compute_fft.jl:30
[5] process_raw!(::SeisData, ::Float64) at /Users/dmikesell/.julia/packages/SeisNoise/I9Knu/src/compute_fft.jl:21
[6] top-level scope at none:0

At the moment, resample! is quite slow on the gpu. Optimizing code therefore requires designing a workflow in which data is resampled before being moved to the gpu so other preprocessing steps can enjoy the gpu speedup. This may not be ideal in all cases. The example below should demonstrate the issue.

# choose original fs and new fs
fs = 100.
resample_fs = 10.

# generate some 'data'
data = rand(500,500)

# resample on cpu
C = CorrData(corr=data,fs=fs)
cpu_time = @elapsed resample!(C,resample_fs)

# resample on gpu
C = CorrData(corr=data,fs=fs) |> gpu
gpu_time = @elapsed resample!(C,resample_fs)

# print results
print("CPU time: ", cpu_time, " s\n")
print("GPU time: ", gpu_time, " s\n")

CPU time: 0.029605658 s
GPU time: 108.995494081 s

Hi Tim,

smooth!() in src/tools.jl might cause unexpected array size change for 1d Array as follows:

using SeisNoise
t = collect(0.0:0.05:60.0)
A = sin.(0.2*pi.*t)

# first time
smooth!(A);
len1 = length(A)
# second time
smooth!(A);
len2 = length(A)
# third time
smooth!(A);
len3 = length(A)

OUTPUT:

julia> println((len1,len2,len3))
(1208, 1215, 1222)

I found that the movingaverage() could cause this issue; the modification of A::AbstractArray seems to affect the returned value (although it shouldn't be without !). So this might be Julia-related issue.

For the moment, I modified them as follows:

function smooth_debug!(A::AbstractArray, half_win::Int=3)
    if ndims(A) == 1
        return movingaverage_debug(A,half_win)
    end

    Nrows, Ncols = size(A)

    for ii = 1:Ncols
        A[:,ii] .= movingaverage_debug(A[:,ii],half_win)
    end
    return nothing
end

function movingaverage_debug(U::AbstractArray, half_win::Int=3)
    #NOTE: avoid to rewrite A
    A = deepcopy(U)
    #----------#
    prepend!(A,A[1:half_win])
    append!(A,A[end-half_win:end])
    N = length(A)
    B = zeros(eltype(A),N)
    window_len = 2 * half_win + 1
    s = sum(A[1:window_len])
    B[half_win+1] = s
    for ii = half_win+2:N-half_win
      s = s - A[ii-half_win] + A[ii+half_win]
      B[ii] = s
    end
    B ./= window_len
    return B[half_win+1:end-half_win-1]
end

And it's solved:

B = sin.(0.2*pi.*t)
# first time
smooth_debug!(B);
len1 = length(B)
# second time
smooth_debug!(B);
len2 = length(B)
# third time
smooth_debug!(B);
len3 = length(B)

OUTPUT:

julia> println((len1,len2,len3))
(1201, 1201, 1201)

Can you manage that until we figure out the call by reference in Julia?

Kurama

Below is an example. The data length is 600 seconds, while the cc_len is 1800 seconds. The last line fails.

using SeisIO, SeisNoise

fs      = 20. # [Hz] resample frequency
cc_len  = 1800 # [s]
cc_step = 450 # [s]

S1 = get_data("FDSN", "TA.I02A..BHZ", src="IRIS", s=string("2007-04-16"), t=600)
S2 = get_data("FDSN", "TA.I02A..BHZ", src="IRIS", s=string("2007-04-16"), t=600)

process_raw!(S1,fs)
process_raw!(S2,fs)

R = RawData.([S1,S2], cc_len, cc_step)

There error is

ArgumentError: invalid index: nothing of type Nothing
to_index(::Nothing) at indices.jl:273
to_index(::Array{Float64,1}, ::Nothing) at indices.jl:250
to_indices at indices.jl:301 [inlined]
to_indices at indices.jl:298 [inlined]
getindex at abstractarray.jl:981 [inlined]
nearest_start_end(::SeisChannel, ::Int64, ::Int64) at slicing.jl:81
RawData(::SeisData, ::Int64, ::Int64) at RawData.jl:77
_broadcast_getindex at broadcast.jl:630 [inlined]
getindex at broadcast.jl:563 [inlined]
copy(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},Type{RawData},Tuple{Array{SeisData,1},Int64,Int64}}) at broadcast.jl:853
materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,Type{RawData},Tuple{Array{SeisData,1},Int64,Int64}}) at broadcast.jl:819
top-level scope at SeisIO_msr_dump.jl:18

It would be better to check the inputs and tell the user they have screwed up and the requested cc_len is longer than the input data. Therefore, there is nothing to index over in terms of time windows. It took me a little while to figure that error out earlier today. I think an assertion should do it.

@assert length(S1.x[1]) >= cc_len

You may need to check for the minimum length input SeisData channels though first.

When the process_raw! function is applied to NodalData, the following error is produced:

MethodError: no method matching process_raw!(::SeisIO.Nodal.NodalData, ::Float64)
Closest candidates are:
  process_raw!(::SeisData, ::Real; ϕshift) at ~/.julia/packages/SeisNoise/E7Pc3/src/compute_fft.jl:20
  process_raw!(::SeisChannel, ::Real; ϕshift) at ~/.julia/packages/SeisNoise/E7Pc3/src/compute_fft.jl:51

Stacktrace:
 [1] top-level scope
   @ In[10]:1

Applying detrend!, taper! and bandpass! to NodalData works fine, so it may be a simple matter of adding NodalData as an acceptable input type for process_raw!, or defining a new function to specifically handle the NodalData.

Hi Tim,
What do you think about returning the filename in save_corr(). That would not change anything other than give the user the filename. We would like that for the AWS tools we are building. Please let us know if you're open to that.
Thanks.

At the moment, the nodalxcorr function outputs a matrix of correlation functions rather than a CorrData object. This means the result of nodalxcorr has no metadata and requires the user to do their own bookkeeping. Ideally, it might output a single NodalCorrData object, which would have a field NodalCorrData.corr containing the matrix of correlation functions for all channel pairs. The NodalCorrDat object would also contain metadata about which station pairs correspond to each correlation function in the matrix.

In case anyone wants to use nodalxcorr and is confused about bookkeeping, the indices of channel pairs for the matrix of correlation functions can be found using the lines below:

Using Combinatorics

# choose which channels to include in the correlation computation
channel_indices = [1,2,3,4]

# get the channel pairs corresponding to each line of the correlation matrix
channel_pair_indices = [j for j in combinations(channel_indices,2)]

6-element Vector{Vector{Int64}}:
 [1, 2]
 [1, 3]
 [1, 4]
 [2, 3]
 [2, 4]
 [3, 4]

The result is a vector in which each element contains the indices of the two channels used to compute the correlation function in the corresponding row of the matrix returned by nodalxcorr.

I am a freshman in JULIA.I want to know the meaning and usage of each parameter in the read_data () function and get_data() function,what can I do?

When computing the errors, should sx2 first be summed over all time windows at each frequency (yielding an array the same length as the frequency array)? Then, instead of dividing s2x2 by a scalar, it is element-wise division of two arrays: the window-summed sx2 and the current window's s2x2.
It would turn

https://github.com/tclements/Noise.jl/blob/c64f6973ccbc2e84cbaa302bf897062fb38eaec2/src/VelocityChange/MWCS.jl#L152-L165

into

# Calculate the slope with a weighted least square linear regression
# forced through the origin
# Pre-sum sx2 over N windows
sx2 = zeros(length(v))
for ii=1:N sx2 .+= sum(w[:,ii] .* v.^2) end
for ii = 1:N
    model = glm(@formula(Y ~0 + X),DataFrame(X=v,Y=phi[:,ii]),Normal(),
                IdentityLink(),wts=w[:,ii])
    # a much faster, unweighted version of this is: m = v \ phi
    dt[ii] = coef(model)[1]

    # Errors
    # σ^2_ϕ
    e = sum((phi[:,ii] .- dt[ii] .* v).^2) / (length(v) - 1)

    s2x2 = (v .* w[:,ii]).^2
    err[ii] = sqrt(e * sum(s2x2 ./ sx2.^2))
end

Dear Professor Clements,
I am a first-year graduate student from China. I still get some wrongs in SeisIO's read_data fuction.
I can't read the data in sac format for disk like yours in the website
SeisNoise.jl/preprocessing.md at master · tclements/SeisNoise.jl (github.com).
My codes are as follows, the program running result is shown in the attached picture.
julia>using Pkg
julia>Pkg.build()
julia>using SeisIO,SeisNoise,Plots
julia>S=read_data("sac","2022.07.31.00.40.31.YD302.BHE.sac")
#"2022.07.31.00.40.31.YD302.BHE.sac" is the name of a sac file

Here is a minimum working example of a bug related to the input SeisData channel and where or not it has a normal PZResp or a multistage response.

using SeisIO, SeisNoise

fs      = 20. # [Hz] resample frequency
cc_len  = 1800 # [s]
cc_step = 450 # [s]

S1 = get_data("FDSN", "TA.I02A..BHZ", src="IRIS", s=string("2007-04-16"), t=6000, msr=true)
S2 = get_data("FDSN", "TA.I02A..BHZ", src="IRIS", s=string("2007-04-16"), t=6000)

process_raw!(S1,fs)
process_raw!(S2,fs)

R = RawData.([S1,S2], cc_len, cc_step)

The results error is related to an inability to convert object types. It would be nice if RawData() could handle either.

MethodError: Cannot convert an object of type MultiStageResp to an object of type PZResp
Closest candidates are:
convert(::Type{T}, !Matched::LazyJSON.PropertyDicts.PropertyDict) where T at /Users/dmikesell/.julia/packages/LazyJSON/nTpJx/src/PropertyDicts.jl:61
convert(::Type{T}, !Matched::LazyJSON.Object) where T at /Users/dmikesell/.julia/packages/LazyJSON/nTpJx/src/AbstractDict.jl:82
convert(::Type{S}, !Matched::T) where {S, T<:(Union{CategoricalArrays.CategoricalString{R}, CategoricalArrays.CategoricalValue{T,R} where T} where R)} at /Users/dmikesell/.julia/packages/CategoricalArrays/dmrjI/src/value.jl:103
...
RawData(::SeisData, ::Int64, ::Int64) at RawData.jl:84
_broadcast_getindex at broadcast.jl:630 [inlined]
getindex at broadcast.jl:563 [inlined]
copy(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},Type{RawData},Tuple{Array{SeisData,1},Int64,Int64}}) at broadcast.jl:853
materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,Type{RawData},Tuple{Array{SeisData,1},Int64,Int64}}) at broadcast.jl:819
top-level scope at SeisIO_msr_dump.jl:13

Hello Timothy,

I'm trying to use arguments in corr_type's correlate function. I am particularly interested in using coherence and PCC. However, these do not work. For example:
C = correlate (FFT [jj], FFT [kk],maxlag, corr_type= "PCC")

I am also interested in using the phase-weighted stack in the next sentence:
stack!(C, allstack = true)
How include the argument?

Thanks for the feedback

At the moment, resample! is about two orders of magnitude slower when applied to NodalData than to CorrData. The two result in very sightly different but essentially identical outputs. The example below should demonstrate the issue.

# choose original fs and new fs
fs = 100.
resample_fs = 10.

# generate some 'data'
sz = 500
data = rand(Float32,sz,sz)

# resample CorrData with same underlying data
C = CorrData(corr=data,fs=fs)
corr_time = @elapsed SeisNoise.resample!(C,resample_fs)

# resample  NodalData with same underlying data
S = SeisNoise.NodalData()
S.data=data
S.fs=ones(sz).*fs
nodal_time = @elapsed SeisNoise.resample!(S,resample_fs)

# print results
print("CorrData time: ", corr_time, " s\n")
print("NodalData time: ", nodal_time, " s\n")

# verify the results are almost identical
print("Total difference: ",sum(C.corr-S.data))

CorrData time: 0.032703704 s
NodalData time: 3.133497156 s
Total difference: 7.137656e-6

At the moment, it seems to be faster to load data into a CorrData object, resample, and then put the data back into a NodalData object.

/home/Fangfan/Pictures/2023-03-02 10-52-13 的屏幕截图.png

Hi Tim, I'm trying to calculate autocorrelations and correlations between components for noise daily records. The script loads the records of all available days of the year and saves the stacked signals in SAC format. However, I observe that the daily autocorrelations do not present any variation between them (the spectrum is identical ). Therefore, I doubt whether the correlations between components are being calculated correctly. I appreciate any feedback about the script I'm using. Attach the Julia script and a day noise record in the three components

ArchivoComprimido.zip



.

https://github.com/tclements/Noise.jl/blob/ca996c929b7bbdcef03ca70a40701919a6193211/src/ArrayFuncs.jl#L20-L22

Should this be either coeff = pinv(A) * X or coeff = inv(A' * A) * A' * X? This may affect detrending functions.