88 lines
3.4 KiB
Go
88 lines
3.4 KiB
Go
|
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
|||
|
|
// SPDX-License-Identifier: MIT
|
|||
|
|
|
|||
|
|
package signal
|
|||
|
|
|
|||
|
|
import (
|
|||
|
|
"sourcedock.dev/petrbalvin/tensor/internal/base"
|
|||
|
|
"sourcedock.dev/petrbalvin/tensor/internal/core"
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
// Zero-phase filtering. A causal IIR filter delays every
|
|||
|
|
// feature by its group delay; running the same filter backwards over
|
|||
|
|
// its own output doubles the magnitude response and cancels the phase
|
|||
|
|
// exactly, because the reversed pass carries the conjugated transfer
|
|||
|
|
// function. The cost is the edge question: the first pass starts from
|
|||
|
|
// rest against a signal that did not, and its start-up transient would
|
|||
|
|
// otherwise bleed into the answer's first samples.
|
|||
|
|
|
|||
|
|
// Filtfilt filters the rank-1 real signal x forwards and backwards
|
|||
|
|
// with the same transfer function the filter designs hand out (b over
|
|||
|
|
// a, a[0] non-zero and normalised away), and returns a result the
|
|||
|
|
// length of x whose magnitude response is the square of the one-pass
|
|||
|
|
// filter's and whose phase is zero: a sinusoid in the passband comes
|
|||
|
|
// out aligned with its input, not lagged, and the group delay at
|
|||
|
|
// every frequency is 0 samples. The dtype contract is FilterApply's:
|
|||
|
|
// float32 and float64 keep their dtype, other real dtypes widen to
|
|||
|
|
// float64.
|
|||
|
|
//
|
|||
|
|
// Edge initialisation: both ends of the signal are extended by an
|
|||
|
|
// even reflection of pad = 3·(nfilt−1) samples (nfilt the longer
|
|||
|
|
// coefficient list), the edge sample not repeated, so the extension
|
|||
|
|
// is continuous in value at both seams. The pad length is the
|
|||
|
|
// standard three times the filter's memory: for a stable filter the
|
|||
|
|
// start-up transient decays like the impulse response tail, which
|
|||
|
|
// 3·(nfilt−1) samples of a direct-form kernel drive below the rounding
|
|||
|
|
// floor for every design this package produces. The forward pass runs
|
|||
|
|
// over the padded signal, the signal is time-reversed, the second
|
|||
|
|
// pass runs, and the pad region of both ends is cropped away, so the
|
|||
|
|
// residual seam error, a slope kink the reflection cannot hide from a
|
|||
|
|
// filter that differentiates, stays outside the returned range. A
|
|||
|
|
// signal no longer than the pad (length ≤ 3·(nfilt−1)) leaves nothing
|
|||
|
|
// to return once both transient regions are excluded and is refused.
|
|||
|
|
func Filtfilt(b, a []float64, x *core.Array) (*core.Array, error) {
|
|||
|
|
const name = "Filtfilt"
|
|||
|
|
bc, ac, out, err := filterPrepare(name, b, a, x)
|
|||
|
|
if err != nil {
|
|||
|
|
return nil, err
|
|||
|
|
}
|
|||
|
|
n := x.Len()
|
|||
|
|
nfilt := max(len(b), len(a))
|
|||
|
|
pad := 3 * (nfilt - 1)
|
|||
|
|
if n <= pad {
|
|||
|
|
return nil, base.Errf("%s: the length-%d signal must exceed the pad length 3·(%d−1) = %d this filter needs",
|
|||
|
|
name, n, nfilt, pad)
|
|||
|
|
}
|
|||
|
|
src := widenFloats(x)
|
|||
|
|
// The reflected extension: left pad holds x[pad], x[pad−1], …,
|
|||
|
|
// x[1] in reverse order, the right pad mirrors it, and the seam at
|
|||
|
|
// either end repeats no sample (whole-sample symmetric even
|
|||
|
|
// reflection).
|
|||
|
|
padded := make([]float64, n+2*pad)
|
|||
|
|
copy(padded[pad:pad+n], src)
|
|||
|
|
for k := range pad {
|
|||
|
|
padded[k] = src[pad-k]
|
|||
|
|
padded[pad+n+k] = src[n-2-k]
|
|||
|
|
}
|
|||
|
|
p := len(padded)
|
|||
|
|
fwd := make([]float64, p)
|
|||
|
|
filterSweep(bc, ac, padded, fwd)
|
|||
|
|
rev := make([]float64, p)
|
|||
|
|
for i := range p {
|
|||
|
|
rev[i] = fwd[p-1-i]
|
|||
|
|
}
|
|||
|
|
back := make([]float64, p)
|
|||
|
|
filterSweep(bc, ac, rev, back)
|
|||
|
|
outF := out.RawFloats()
|
|||
|
|
if outF == nil || out.Strided() {
|
|||
|
|
for t := range n {
|
|||
|
|
out.SetFloatAt(t, back[p-1-pad-t])
|
|||
|
|
}
|
|||
|
|
} else {
|
|||
|
|
for t := range n {
|
|||
|
|
outF[t] = back[p-1-pad-t]
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return out, nil
|
|||
|
|
}
|