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
|
||
}
|