From d7cb166f3d8671713a33d17b1b4b459214e10501 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 31 Jan 2024 19:49:00 -0500 Subject: [PATCH] dsp: fix filt to filtfilt --- Project.toml | 2 +- .../Preconditioners/DataPreconditioners.jl | 14 +++++++++++--- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index df06dc4c2..696221da6 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ JOLI = "0.7, 0.8" OrderedCollections = "1.5" PyCall = "1.18, 1.90, 1.91, 1.62" Requires = "1" -SegyIO = "0.7.7 - 0.8.3" +SegyIO = "0.7.7 - 0.8.3, >= 0.8.5" TimerOutputs = "0.5" julia = "1.6" diff --git a/src/TimeModeling/Preconditioners/DataPreconditioners.jl b/src/TimeModeling/Preconditioners/DataPreconditioners.jl index fa9608e6b..fe5d5dd13 100644 --- a/src/TimeModeling/Preconditioners/DataPreconditioners.jl +++ b/src/TimeModeling/Preconditioners/DataPreconditioners.jl @@ -175,6 +175,14 @@ conj(I::FrequencyFilter{T}) where T = I adjoint(I::FrequencyFilter{T}) where T = I transpose(I::FrequencyFilter{T}) where T = I +function tracefilt!(x, y, ypad, filter) + n = length(y) + ypad[1:n-1] .= view(y, n:-1:2) + ypad[n:end] .= y + x .= filtfilt(filter, ypad)[n:end] + nothing +end + function filter!(dout::AbstractArray{T, N}, din::AbstractArray{T, N}, dt::T; fmin=T(0.01), fmax=T(100)) where {T, N} if fmin == 0 responsetype = Lowpass(fmax; fs=1e3/dt) @@ -183,9 +191,9 @@ function filter!(dout::AbstractArray{T, N}, din::AbstractArray{T, N}, dt::T; fmi else responsetype = Bandpass(fmin, fmax; fs=1e3/dt) end - designmethod = Butterworth(5) - tracefilt!(x, y) = filt!(x, digitalfilter(responsetype, designmethod), y) - map(i-> tracefilt!(selectdim(dout, N, i), selectdim(din, N, i)), 1:size(dout, 2)) + filter = digitalfilter(responsetype, Butterworth(5)) + ytmp = zeros(T, size(din, 1)*2-1) + map(i-> tracefilt!(selectdim(dout, N, i), selectdim(din, N, i), ytmp, filter), 1:size(dout, 2)) end filter_data(Din::judiVector; fmin=0, fmax=Inf) = judiFilter(Din.geometry, fmin, fmax)*Din