From 002a315d414e5ac2ba264e939336f8b2b1770563 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Sat, 21 Oct 2023 09:16:13 +0200
Subject: [PATCH] add method to `DSP.filt` and `filtfilt`

---
 lib/ControlSystemsBase/src/dsp.jl       | 28 ++++++++++++++++++++++++-
 lib/ControlSystemsBase/test/test_dsp.jl | 13 ++++++++++++
 2 files changed, 40 insertions(+), 1 deletion(-)

diff --git a/lib/ControlSystemsBase/src/dsp.jl b/lib/ControlSystemsBase/src/dsp.jl
index c580926bc..c8c5cde44 100644
--- a/lib/ControlSystemsBase/src/dsp.jl
+++ b/lib/ControlSystemsBase/src/dsp.jl
@@ -36,4 +36,30 @@ end
 function zpk(p::DSP.ZeroPoleGain, h::Real)
     z,p,k = p.z, p.p, p.k
     zpk(z, p, k, h)
-end
\ No newline at end of file
+end
+
+"""
+    DSP.filt(P::ControlSystemsBase.TransferFunction, u)
+
+Use a transfer function `P` to filter a signal `u`. This is equivalent to `lsim(P, u').y[:]`, but may be more efficient for single-input, single-output systems when the state sequence is not needed.
+"""
+function DSP.filt(P::ControlSystemsBase.TransferFunction, u, args...)
+    issiso(P) || error("Only single-input, single-output systems are supported in filt, call lsim instead.")
+    b, a = numvec(P)[], denvec(P)[]
+    nb, na = length(b), length(a)
+    if nb <= na
+        b = [zeros(na - nb); b]
+    end
+    DSP.filt(b, a, u, args...)
+end
+
+
+function DSP.filtfilt(P::ControlSystemsBase.TransferFunction, u, args...)
+    issiso(P) || error("Only single-input, single-output systems are supported in filtfilt.")
+    b, a = numvec(P)[], denvec(P)[]
+    nb, na = length(b), length(a)
+    if nb <= na
+        b = [zeros(na - nb); b]
+    end
+    DSP.filtfilt(b, a, u, args...)
+end
diff --git a/lib/ControlSystemsBase/test/test_dsp.jl b/lib/ControlSystemsBase/test/test_dsp.jl
index 532a460bb..4bbf9e556 100644
--- a/lib/ControlSystemsBase/test/test_dsp.jl
+++ b/lib/ControlSystemsBase/test/test_dsp.jl
@@ -26,5 +26,18 @@
     uf = DSP.filt(f, u)
     uls = lsim(fcs, u').y'
     @test uf ≈ uls
+
+
+
+
+
+
+    u = randn(100)
+    res = lsim(Gd, u')
+    yf = res.y
+    @test yf' ≈ DSP.filt(Gd, u)
+
+    # The filtfilt function has a complicated initialization of the filter state, making it hard to compare to twice lsim
+    @test length(DSP.filtfilt(Gd, u)) == length(u)
     
 end