diff --git a/examples/iir.cpp b/examples/iir.cpp
index 8ece8295..2c5fbeee 100644
--- a/examples/iir.cpp
+++ b/examples/iir.cpp
@@ -14,91 +14,84 @@ int main()
 {
     println(library_version());
 
-    constexpr size_t maxorder = 32;
-
     const std::string options = "phaseresp=True, log_freq=True, freq_dB_lim=(-160, 10), padwidth=8192";
 
     univector<fbase, 1024> output;
     {
-        zpk<fbase> filt                       = iir_lowpass(bessel<fbase>(24), 1000, 48000);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt = iir_lowpass(bessel<fbase>(24), 1000, 48000);
+        output          = iir(unitimpulse(), filt);
     }
     plot_save("bessel_lowpass24", output, options + ", title='24th-order Bessel filter, lowpass 1khz'");
 
     {
-        zpk<fbase> filt                       = iir_lowpass(bessel<fbase>(12), 1000, 48000);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt = iir_lowpass(bessel<fbase>(12), 1000, 48000);
+        output          = iir(unitimpulse(), filt);
     }
     plot_save("bessel_lowpass12", output, options + ", title='12th-order Bessel filter, lowpass 1khz'");
 
     {
-        zpk<fbase> filt                       = iir_lowpass(bessel<fbase>(6), 1000, 48000);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt = iir_lowpass(bessel<fbase>(6), 1000, 48000);
+        output          = iir(unitimpulse(), filt);
     }
     plot_save("bessel_lowpass6", output, options + ", title='6th-order Bessel filter, lowpass 1khz'");
 
     {
-        zpk<fbase> filt                       = iir_lowpass(butterworth<fbase>(24), 1000, 48000);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt = iir_lowpass(butterworth<fbase>(24), 1000, 48000);
+        output          = iir(unitimpulse(), filt);
     }
     plot_save("butterworth_lowpass24", output,
               options + ", title='24th-order Butterworth filter, lowpass 1khz'");
 
     {
-        zpk<fbase> filt                       = iir_lowpass(butterworth<fbase>(12), 1000, 48000);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt = iir_lowpass(butterworth<fbase>(12), 1000, 48000);
+        output          = iir(unitimpulse(), filt);
     }
     plot_save("butterworth_lowpass12", output,
               options + ", title='12th-order Butterworth filter, lowpass 1khz'");
 
     {
-        zpk<fbase> filt                       = iir_highpass(butterworth<fbase>(12), 1000, 48000);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt       = iir_highpass(butterworth<fbase>(12), 1000, 48000);
+        iir_params<fbase> bqs = to_sos(filt); // to_sos is expensive, keep iir_params if reused
+        output                = iir(unitimpulse(), bqs);
     }
     plot_save("butterworth_highpass12", output,
               options + ", title='12th-order Butterworth filter, highpass 1khz'");
 
     {
-        zpk<fbase> filt                       = iir_bandpass(butterworth<fbase>(12), 0.1, 0.2);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt       = iir_bandpass(butterworth<fbase>(12), 0.1, 0.2);
+        iir_params<fbase> bqs = to_sos(filt); // to_sos is expensive, keep iir_params if reused
+        output                = iir(unitimpulse(), bqs);
     }
     plot_save("butterworth_bandpass12", output,
               options + ", title='12th-order Butterworth filter, bandpass'");
 
     {
-        zpk<fbase> filt                       = iir_bandstop(butterworth<fbase>(12), 0.1, 0.2);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt       = iir_bandstop(butterworth<fbase>(12), 0.1, 0.2);
+        iir_params<fbase> bqs = to_sos(filt); // to_sos is expensive, keep iir_params if reused
+        output                = iir(unitimpulse(), bqs);
     }
     plot_save("butterworth_bandstop12", output,
               options + ", title='12th-order Butterworth filter, bandstop'");
 
     {
-        zpk<fbase> filt                       = iir_bandpass(butterworth<fbase>(4), 0.005, 0.9);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt       = iir_bandpass(butterworth<fbase>(4), 0.005, 0.9);
+        iir_params<fbase> bqs = to_sos(filt); // to_sos is expensive, keep iir_params if reused
+        output                = iir(unitimpulse(), bqs);
     }
     plot_save("butterworth_bandpass4", output, options + ", title='4th-order Butterworth filter, bandpass'");
 
     {
-        zpk<fbase> filt                       = iir_lowpass(chebyshev1<fbase>(8, 2), 0.09);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt       = iir_lowpass(chebyshev1<fbase>(8, 2), 0.09);
+        iir_params<fbase> bqs = to_sos(filt); // to_sos is expensive, keep iir_params if reused
+        output                = iir(unitimpulse(), bqs);
     }
     plot_save("chebyshev1_lowpass8", output,
               options + ", title='8th-order Chebyshev type I filter, lowpass'");
 
     {
-        zpk<fbase> filt                       = iir_lowpass(chebyshev2<fbase>(8, 80), 0.09);
-        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
-        output                                = biquad<maxorder>(bqs, unitimpulse());
+        zpk<fbase> filt       = iir_lowpass(chebyshev2<fbase>(8, 80), 0.09);
+        iir_params<fbase> bqs = to_sos(filt); // to_sos is expensive, keep iir_params if reused
+        output                = iir(unitimpulse(), filt);
     }
     plot_save("chebyshev2_lowpass8", output,
               options + ", title='8th-order Chebyshev type II filter, lowpass'");
diff --git a/include/kfr/base/state_holder.hpp b/include/kfr/base/state_holder.hpp
index 763a6445..e2d6a05c 100644
--- a/include/kfr/base/state_holder.hpp
+++ b/include/kfr/base/state_holder.hpp
@@ -15,13 +15,20 @@ namespace kfr
 {
 
 template <typename T, bool Stateless>
-struct state_holder
+struct state_holder;
+
+template <typename T>
+struct state_holder<T, false>
 {
+    static_assert(!std::is_const_v<T>, "state_holder: T must not be const");
+
     constexpr state_holder()                    = delete;
     constexpr state_holder(const state_holder&) = default;
     constexpr state_holder(state_holder&&)      = default;
-    constexpr state_holder(const T& state) CMT_NOEXCEPT : s(state) {}
-    constexpr state_holder(std::reference_wrapper<const T> state) CMT_NOEXCEPT : s(state) {}
+    constexpr state_holder(T state) CMT_NOEXCEPT : s(std::move(state)) {}
+    constexpr state_holder(std::reference_wrapper<T> state)       = delete;
+    constexpr state_holder(std::reference_wrapper<const T> state) = delete;
+    constexpr state_holder(state_holder<T, true> stateless) : s(*stateless) {}
     T s;
 
     const T* operator->() const { return &s; }
@@ -33,17 +40,28 @@ struct state_holder
 template <typename T>
 struct state_holder<T, true>
 {
-    constexpr state_holder()                    = delete;
-    constexpr state_holder(const state_holder&) = default;
-    constexpr state_holder(state_holder&&)      = default;
-    constexpr state_holder(T& state) CMT_NOEXCEPT : s(state) {}
-    constexpr state_holder(std::reference_wrapper<T> state) CMT_NOEXCEPT : s(state) {}
-    T& s;
+    static_assert(!std::is_const_v<T>, "state_holder: T must not be const");
+    
+    constexpr state_holder()                            = delete;
+    constexpr state_holder(const state_holder&)         = default;
+    constexpr state_holder(state_holder&&)              = default;
+    constexpr state_holder(T state) CMT_NOEXCEPT        = delete;
+    constexpr state_holder(const T& state) CMT_NOEXCEPT = delete;
+    constexpr state_holder(T& state) CMT_NOEXCEPT       = delete;
+    constexpr state_holder(T&& state) CMT_NOEXCEPT      = delete;
+    constexpr state_holder(std::reference_wrapper<T> state) CMT_NOEXCEPT : s(&state.get()) {}
+    T* s;
 
-    const T* operator->() const { return &s; }
-    T* operator->() { return &s; }
-    const T& operator*() const { return s; }
-    T& operator*() { return s; }
+    const T* operator->() const { return s; }
+    T* operator->() { return s; }
+    const T& operator*() const { return *s; }
+    T& operator*() { return *s; }
 };
 
+static_assert(std::is_copy_constructible_v<state_holder<float, true>>);
+static_assert(std::is_copy_constructible_v<state_holder<float, false>>);
+
+static_assert(std::is_move_constructible_v<state_holder<float, true>>);
+static_assert(std::is_move_constructible_v<state_holder<float, false>>);
+
 } // namespace kfr
diff --git a/include/kfr/cometa.hpp b/include/kfr/cometa.hpp
index bb17b5eb..70b7b060 100644
--- a/include/kfr/cometa.hpp
+++ b/include/kfr/cometa.hpp
@@ -1433,6 +1433,18 @@ struct has_data_size_impl<T,
 {
 };
 
+template <typename, typename = void>
+struct has_data_size_free_impl : std::false_type
+{
+};
+
+template <typename T>
+struct has_data_size_free_impl<
+    T, std::void_t<decltype(std::size(std::declval<T>())), decltype(std::data(std::declval<T>()))>>
+    : std::true_type
+{
+};
+
 template <typename, typename Fallback, typename = void>
 struct value_type_impl
 {
@@ -1453,8 +1465,7 @@ template <typename T>
 constexpr inline bool has_data_size = details::has_data_size_impl<std::decay_t<T>>::value;
 
 #define CMT_HAS_DATA_SIZE(CONTAINER)                                                                         \
-    std::void_t<decltype(std::size(std::declval<CONTAINER>())),                                              \
-                decltype(std::data(std::declval<CONTAINER>()))>* = nullptr
+    std::enable_if_t<details::has_data_size_free_impl<CONTAINER>::value>* = nullptr
 
 template <typename T>
 using value_type_of = typename std::decay_t<T>::value_type;
diff --git a/include/kfr/dsp/biquad.hpp b/include/kfr/dsp/biquad.hpp
index d574eac1..ef26e887 100644
--- a/include/kfr/dsp/biquad.hpp
+++ b/include/kfr/dsp/biquad.hpp
@@ -29,12 +29,21 @@
 #include "../base/handle.hpp"
 #include "../simd/impl/function.hpp"
 #include "../simd/operators.hpp"
+#include "../base/state_holder.hpp"
 #include "../simd/vec.hpp"
 #include "../testo/assert.hpp"
 
 namespace kfr
 {
 
+constexpr inline size_t maximum_iir_order    = 128;
+constexpr inline size_t maximum_biquad_count = maximum_iir_order / 2;
+
+namespace internal_generic
+{
+constexpr inline auto biquad_sizes = csize<1> << csizeseq<ilog2(maximum_biquad_count) + 1>;
+}
+
 enum class biquad_type
 {
     lowpass,
@@ -51,26 +60,26 @@ enum class biquad_type
  * @brief Structure for holding biquad filter coefficients.
  */
 template <typename T>
-struct biquad_params
+struct biquad_section
 {
     template <typename U>
-    constexpr biquad_params(const biquad_params<U>& bq) CMT_NOEXCEPT : a0(static_cast<T>(bq.a0)),
-                                                                       a1(static_cast<T>(bq.a1)),
-                                                                       a2(static_cast<T>(bq.a2)),
-                                                                       b0(static_cast<T>(bq.b0)),
-                                                                       b1(static_cast<T>(bq.b1)),
-                                                                       b2(static_cast<T>(bq.b2))
+    constexpr biquad_section(const biquad_section<U>& bq) CMT_NOEXCEPT : a0(static_cast<T>(bq.a0)),
+                                                                         a1(static_cast<T>(bq.a1)),
+                                                                         a2(static_cast<T>(bq.a2)),
+                                                                         b0(static_cast<T>(bq.b0)),
+                                                                         b1(static_cast<T>(bq.b1)),
+                                                                         b2(static_cast<T>(bq.b2))
     {
     }
 
     static_assert(std::is_floating_point_v<T>, "T must be a floating point type");
-    constexpr biquad_params() CMT_NOEXCEPT : a0(1), a1(0), a2(0), b0(1), b1(0), b2(0) {}
-    constexpr biquad_params(T a0, T a1, T a2, T b0, T b1, T b2) CMT_NOEXCEPT : a0(a0),
-                                                                               a1(a1),
-                                                                               a2(a2),
-                                                                               b0(b0),
-                                                                               b1(b1),
-                                                                               b2(b2)
+    constexpr biquad_section() CMT_NOEXCEPT : a0(1), a1(0), a2(0), b0(1), b1(0), b2(0) {}
+    constexpr biquad_section(T a0, T a1, T a2, T b0, T b1, T b2) CMT_NOEXCEPT : a0(a0),
+                                                                                a1(a1),
+                                                                                a2(a2),
+                                                                                b0(b0),
+                                                                                b1(b1),
+                                                                                b2(b2)
     {
     }
     T a0;
@@ -79,14 +88,14 @@ struct biquad_params
     T b0;
     T b1;
     T b2;
-    biquad_params<T> normalized_a0() const
+    biquad_section<T> normalized_a0() const
     {
         vec<T, 5> v{ a1, a2, b0, b1, b2 };
         v = v / a0;
         return { T(1.0), v[0], v[1], v[2], v[3], v[4] };
     }
-    biquad_params<T> normalized_b0() const { return { a0, a1, a2, T(1.0), b1 / b0, b2 / b0 }; }
-    biquad_params<T> normalized_all() const { return normalized_a0().normalized_b0(); }
+    biquad_section<T> normalized_b0() const { return { a0, a1, a2, T(1.0), b1 / b0, b2 / b0 }; }
+    biquad_section<T> normalized_all() const { return normalized_a0().normalized_b0(); }
 };
 
 template <typename T, size_t filters>
@@ -98,8 +107,8 @@ struct biquad_state
     constexpr biquad_state() CMT_NOEXCEPT : s1(0), s2(0), out(0) {}
 };
 
-template <typename T, size_t filters>
-struct biquad_block
+template <typename T, size_t filters = tag_dynamic_vector>
+struct iir_params
 {
     vec<T, filters> a1;
     vec<T, filters> a2;
@@ -107,10 +116,11 @@ struct biquad_block
     vec<T, filters> b1;
     vec<T, filters> b2;
 
-    constexpr biquad_block() CMT_NOEXCEPT : a1(0), a2(0), b0(1), b1(0), b2(0) {}
-    CMT_GNU_CONSTEXPR biquad_block(const biquad_params<T>* bq, size_t count) CMT_NOEXCEPT
+    constexpr iir_params() CMT_NOEXCEPT : a1(0), a2(0), b0(1), b1(0), b2(0) {}
+    CMT_GNU_CONSTEXPR iir_params(const biquad_section<T>* bq, size_t count)
     {
-        count = count > filters ? filters : count;
+        KFR_LOGIC_CHECK(count <= filters, "iir_params: too many biquad sections");
+        count = const_min(filters, count);
         for (size_t i = 0; i < count; i++)
         {
             a1[i] = bq[i].a1;
@@ -129,49 +139,123 @@ struct biquad_block
         }
     }
 
-    template <size_t count>
-    constexpr biquad_block(const biquad_params<T> (&bq)[count]) CMT_NOEXCEPT : biquad_block(bq, count)
+    CMT_GNU_CONSTEXPR iir_params(const biquad_section<T>& one) CMT_NOEXCEPT : iir_params(&one, 1) {}
+
+    template <typename Container, CMT_HAS_DATA_SIZE(Container)>
+    constexpr iir_params(Container&& cont) CMT_NOEXCEPT : iir_params(std::data(cont), std::size(cont))
+    {
+    }
+};
+
+template <typename T>
+struct iir_params<T, tag_dynamic_vector> : public std::vector<biquad_section<T>>
+{
+    using base = std::vector<biquad_section<T>>;
+
+    iir_params()                  = default;
+    iir_params(const iir_params&) = default;
+    iir_params(iir_params&&)      = default;
+
+    iir_params(size_t count) : base(count) {}
+
+    iir_params(const biquad_section<T>* bq, size_t count) CMT_NOEXCEPT : base(bq, bq + count) {}
+
+    iir_params(const biquad_section<T>& one) CMT_NOEXCEPT : iir_params(&one, 1) {}
+
+    iir_params(std::vector<biquad_section<T>>&& sections) CMT_NOEXCEPT : base(std::move(sections)) {}
+
+    template <typename Container, CMT_HAS_DATA_SIZE(Container)>
+    constexpr iir_params(Container&& cont) CMT_NOEXCEPT : iir_params(std::data(cont), std::size(cont))
+    {
+    }
+
+    template <size_t filters>
+    iir_params(const iir_params<T, filters>& params) : base(filters)
     {
-        static_assert(count <= filters, "count > filters");
+        for (size_t i = 0; i < filters; ++i)
+        {
+            this->operator[](i).a0 = T(1);
+            this->operator[](i).a1 = params.a1[i];
+            this->operator[](i).a2 = params.a2[i];
+            this->operator[](i).b0 = params.b0[i];
+            this->operator[](i).b1 = params.b1[i];
+            this->operator[](i).b2 = params.b2[i];
+        }
     }
 };
 
+template <typename T, size_t Size>
+iir_params(const std::array<T, Size>&) -> iir_params<T, Size>;
+template <typename T, size_t Size>
+iir_params(const univector<T, Size>) -> iir_params<T, Size>;
+template <typename T, size_t Size>
+iir_params(const biquad_section<T> (&)[Size]) -> iir_params<T, Size>;
+template <typename T>
+iir_params(const biquad_section<T>&) -> iir_params<T, 1>;
+template <typename T>
+iir_params(const std::vector<biquad_section<T>>&) -> iir_params<T, tag_dynamic_vector>;
+template <typename T>
+iir_params(std::vector<biquad_section<T>>&&) -> iir_params<T, tag_dynamic_vector>;
+
+template <typename T, size_t filters>
+struct iir_state
+{
+    static_assert(filters >= 1 && filters <= maximum_biquad_count, "Incorrect number of biquad filters");
+
+    iir_params<T, filters> params;
+
+    template <typename... Args,
+              std::enable_if_t<std::is_constructible_v<iir_params<T, filters>, Args...>>* = nullptr>
+    iir_state(Args&&... args) : params(std::forward<Args>(args)...)
+    {
+    }
+
+    biquad_state<T, filters> state;
+    biquad_state<T, filters> saved_state;
+    size_t block_end = 0;
+};
+
+template <typename T, size_t filters>
+iir_state(const iir_params<T, filters>&) -> iir_state<T, filters>;
+template <typename T, size_t filters>
+iir_state(iir_params<T, filters>&&) -> iir_state<T, filters>;
+
 inline namespace CMT_ARCH_NAME
 {
 
-template <size_t filters, typename T, typename E1>
-struct expression_biquads_l : public expression_with_traits<E1>
+template <size_t filters, typename T, typename E1, bool Stateless = false>
+struct expression_iir_l : public expression_with_traits<E1>
 {
     using value_type = T;
 
-    expression_biquads_l(const biquad_block<T, filters>& bq, E1&& e1)
-        : expression_with_traits<E1>(std::forward<E1>(e1)), bq(bq)
+    expression_iir_l(E1&& e1, state_holder<iir_state<T, filters>, Stateless> state)
+        : expression_with_traits<E1>(std::forward<E1>(e1)), state(std::move(state))
     {
     }
-    biquad_block<T, filters> bq;
-    mutable biquad_state<T, filters> state;
+
+    mutable state_holder<iir_state<T, filters>, Stateless> state;
 };
 
-template <size_t filters, typename T, typename E1>
-struct expression_biquads : expression_with_traits<E1>
+template <size_t filters, typename T, typename E1, bool Stateless = false>
+struct expression_iir : expression_with_traits<E1>
 {
     using value_type = T;
 
-    expression_biquads(const biquad_block<T, filters>& bq, E1&& e1)
-        : expression_with_traits<E1>(std::forward<E1>(e1)), bq(bq), block_end(0)
+    expression_iir(E1&& e1, state_holder<iir_state<T, filters>, Stateless> state)
+        : expression_with_traits<E1>(std::forward<E1>(e1)), state(std::move(state))
     {
     }
 
-    biquad_block<T, filters> bq;
-
-    mutable biquad_state<T, filters> state;
-    mutable biquad_state<T, filters> saved_state;
-    mutable size_t block_end;
+    mutable state_holder<iir_state<T, filters>, Stateless> state;
 };
 
+namespace internal
+{
+
 template <size_t filters, typename T>
-KFR_INTRINSIC T biquad_process(vec<T, filters>& out, const biquad_block<T, filters>& bq,
-                               biquad_state<T, filters>& state, T in0, const vec<T, filters>& delayline)
+KFR_INTRINSIC T biquad_process(vec<T, filters>& out, const iir_params<T, filters>& bq,
+                               biquad_state<T, filters>& state, identity<T> in0,
+                               const vec<T, filters>& delayline)
 {
     vec<T, filters> in = insertleft(in0, delayline);
     out                = bq.b0 * in + state.s1;
@@ -179,94 +263,168 @@ KFR_INTRINSIC T biquad_process(vec<T, filters>& out, const biquad_block<T, filte
     state.s2           = bq.b2 * in - bq.a2 * out;
     return out[filters - 1];
 }
-
-template <size_t filters, typename T, typename E1, size_t N>
-KFR_INTRINSIC vec<T, N> get_elements(const expression_biquads_l<filters, T, E1>& self, shape<1> index,
-                                     axis_params<0, N> t)
+template <size_t filters, typename T, size_t N>
+KFR_INTRINSIC vec<T, N> biquad_process(iir_state<T, filters>& state, const vec<T, N>& in,
+                                       size_t save_state_after = static_cast<size_t>(-1))
 {
-    const vec<T, N> in = get_elements(self.first(), index, t);
     vec<T, N> out;
-
-    CMT_LOOP_UNROLL
-    for (size_t i = 0; i < N; i++)
+    if (CMT_LIKELY(save_state_after == static_cast<size_t>(-1)))
     {
-        out[i] = biquad_process(self.state.out, self.bq, self.state, in[i], self.state.out);
+        CMT_LOOP_UNROLL
+        for (size_t i = 0; i < N; i++)
+        {
+            out[i] = biquad_process(state.state.out, state.params, state.state, in[i], state.state.out);
+        }
+    }
+    else
+    {
+        for (size_t i = 0; i < save_state_after; i++)
+        {
+            out[i] = biquad_process(state.state.out, state.params, state.state, in[i], state.state.out);
+        }
+        state.saved_state = state.state;
+        for (size_t i = save_state_after; i < N; i++)
+        {
+            out[i] = biquad_process(state.state.out, state.params, state.state, in[i], state.state.out);
+        }
     }
-
     return out;
 }
+} // namespace internal
+
+template <size_t filters, typename T, typename E1, size_t N>
+KFR_INTRINSIC vec<T, N> get_elements(const expression_iir_l<filters, T, E1>& self, shape<1> index,
+                                     axis_params<0, N> t)
+{
+    const vec<T, N> in = get_elements(self.first(), index, t);
+    return internal::biquad_process(*self.state, in);
+}
 
+template <typename T, typename E1>
+KFR_INTRINSIC void begin_pass(const expression_iir<1, T, E1>&, shape<1>, shape<1>)
+{
+}
 template <size_t filters, typename T, typename E1>
-KFR_INTRINSIC void begin_pass(const expression_biquads<filters, T, E1>& self, shape<1> start, shape<1> stop)
+KFR_INTRINSIC void begin_pass(const expression_iir<filters, T, E1>& self, shape<1> start, shape<1> stop)
 {
-    size_t size    = stop.front();
-    self.block_end = size;
+    size_t size           = stop.front();
+    self.state->block_end = size;
+    vec<T, filters - 1> in;
     for (index_t i = 0; i < filters - 1; i++)
     {
-        const vec<T, 1> in = i < size ? get_elements(self.first(), shape<1>{ i }, axis_params_v<0, 1>) : 0;
-        biquad_process(self.state.out, self.bq, self.state, in[0], self.state.out);
+        in[i] = i < size ? get_elements(self.first(), shape<1>{ i }, axis_params_v<0, 1>).front() : 0;
     }
+    internal::biquad_process(*self.state, in);
+}
+
+template <typename T, typename E1>
+KFR_INTRINSIC void end_pass(const expression_iir<1, T, E1>&, shape<1>, shape<1>)
+{
 }
 template <size_t filters, typename T, typename E1>
-KFR_INTRINSIC void end_pass(const expression_biquads<filters, T, E1>& self, shape<1> start, shape<1> stop)
+KFR_INTRINSIC void end_pass(const expression_iir<filters, T, E1>& self, shape<1> start, shape<1> stop)
 {
-    self.state = self.saved_state;
+    self.state->state = self.state->saved_state;
+}
+
+template <typename T, typename E1, size_t N>
+KFR_INTRINSIC vec<T, N> get_elements(const expression_iir<1, T, E1>& self, shape<1> index,
+                                     axis_params<0, N> t)
+{
+    const vec<T, N> in = get_elements(self.first(), index, t);
+    return internal::biquad_process(*self.state, in);
 }
 
 template <size_t filters, typename T, typename E1, size_t N>
-KFR_INTRINSIC vec<T, N> get_elements(const expression_biquads<filters, T, E1>& self, shape<1> index,
+KFR_INTRINSIC vec<T, N> get_elements(const expression_iir<filters, T, E1>& self, shape<1> index,
                                      axis_params<0, N> t)
 {
+    using internal::biquad_process;
     index.front() += filters - 1;
     vec<T, N> out{};
-    if (index.front() + N <= self.block_end)
+    if (index.front() + N <= self.state->block_end)
     {
         const vec<T, N> in = get_elements(self.first(), shape<1>{ index.front() }, t);
 
-        CMT_LOOP_UNROLL
-        for (size_t i = 0; i < N; i++)
-        {
-            out[i] = biquad_process(self.state.out, self.bq, self.state, in[i], self.state.out);
-        }
-        if (index.front() + N == self.block_end)
-            self.saved_state = self.state;
+        out = biquad_process(*self.state, in);
+        if (index.front() + N == self.state->block_end)
+            self.state->saved_state = self.state->state;
     }
-    else if (index.front() >= self.block_end)
+    else if (index.front() >= self.state->block_end)
     {
-        CMT_LOOP_UNROLL
-        for (size_t i = 0; i < N; i++)
-        {
-            out[i] = biquad_process(self.state.out, self.bq, self.state, T(0), self.state.out);
-        }
+        out = biquad_process(*self.state, vec<T, N>(0));
     }
     else
     {
-        size_t i = 0;
-        for (; i < std::min(N, self.block_end - static_cast<size_t>(index.front())); i++)
-        {
-            const vec<T, 1> in =
-                get_elements(self.first(), index.add_at(i, cval<index_t, 0>), axis_params_v<0, 1>);
-            out[i] = biquad_process(self.state.out, self.bq, self.state, in[0], self.state.out);
-        }
-        self.saved_state = self.state;
-        for (; i < N; i++)
-        {
-            out[i] = biquad_process(self.state.out, self.bq, self.state, T(0), self.state.out);
-        }
+        size_t save_at = std::min(N, self.state->block_end - static_cast<size_t>(index.front()));
+        vec<T, N> in;
+        for (size_t i = 0; i < save_at; ++i)
+            in[i] =
+                get_elements(self.first(), index.add_at(i, cval<index_t, 0>), axis_params_v<0, 1>).front();
+        for (size_t i = save_at; i < N; ++i)
+            in[i] = 0;
+        out = biquad_process(*self.state, in, save_at);
     }
     return out;
 }
 
+/**
+ * @brief Returns template expressions that applies biquad filter to the input.
+ * @param e1 Input expression
+ * @param params Biquad coefficients
+ */
+template <size_t filters, typename T, typename E1>
+KFR_FUNCTION expression_iir<filters, T, E1> iir(E1&& e1, iir_params<T, filters> params)
+{
+    return expression_iir<filters, T, E1>(std::forward<E1>(e1), iir_state{ std::move(params) });
+}
+
+/**
+ * @brief Returns template expressions that applies biquad filter to the input.
+ * @param e1 Input expression
+ * @param params Biquad coefficients
+ */
+template <typename T, typename E1>
+KFR_FUNCTION expression_handle<T, 1> iir(E1&& e1, const iir_params<T, tag_dynamic_vector>& params)
+{
+    KFR_LOGIC_CHECK(next_poweroftwo(params.size()) <= maximum_biquad_count, "iir: too many biquad sections");
+    return cswitch(
+        internal_generic::biquad_sizes, next_poweroftwo(params.size()),
+        [&](auto x)
+        {
+            constexpr size_t filters = x;
+            return to_handle(expression_iir<filters, T, E1>(
+                std::forward<E1>(e1), iir_state{ iir_params<T, filters>(params.data(), params.size()) }));
+        },
+        [&] { return to_handle(fixshape(zeros<T>(), fixed_shape<infinite_size>)); });
+}
+
+/**
+ * @brief Returns template expressions that applies biquad filter to the input.
+ * @param bq Biquad coefficients
+ * @param e1 Input expression
+ */
+template <size_t filters, typename T, typename E1>
+KFR_FUNCTION expression_iir<filters, T, E1, true> iir(E1&& e1,
+                                                      std::reference_wrapper<iir_state<T, filters>> state)
+{
+    return expression_iir<filters, T, E1, true>(std::forward<E1>(e1), state);
+}
+
+#define KFR_BIQUAD_DEPRECATED                                                                                \
+    [[deprecated("biquad(param, expr) prototype is deprecated. Use iir(expr, param) with swapped "           \
+                 "arguments")]]
+
 /**
  * @brief Returns template expressions that applies biquad filter to the input.
  * @param bq Biquad coefficients
  * @param e1 Input expression
  */
 template <typename T, typename E1>
-KFR_FUNCTION expression_biquads<1, T, E1> biquad(const biquad_params<T>& bq, E1&& e1)
+KFR_BIQUAD_DEPRECATED KFR_FUNCTION expression_iir<1, T, E1> biquad(const biquad_section<T>& bq, E1&& e1)
 {
-    const biquad_params<T> bqs[1] = { bq };
-    return expression_biquads<1, T, E1>(bqs, std::forward<E1>(e1));
+    const biquad_section<T> bqs[1] = { bq };
+    return expression_iir<1, T, E1>(std::forward<E1>(e1), iir_state{ iir_params{ bqs } });
 }
 
 /**
@@ -276,9 +434,10 @@ KFR_FUNCTION expression_biquads<1, T, E1> biquad(const biquad_params<T>& bq, E1&
  * @note This implementation introduces delay of N - 1 samples, where N is the filter count.
  */
 template <size_t filters, typename T, typename E1>
-KFR_FUNCTION expression_biquads_l<filters, T, E1> biquad_l(const biquad_params<T> (&bq)[filters], E1&& e1)
+KFR_BIQUAD_DEPRECATED KFR_FUNCTION expression_iir_l<filters, T, E1> biquad_l(
+    const biquad_section<T> (&bq)[filters], E1&& e1)
 {
-    return expression_biquads_l<filters, T, E1>(bq, std::forward<E1>(e1));
+    return expression_iir_l<filters, T, E1>(std::forward<E1>(e1), iir_state{ iir_params{ bq } });
 }
 
 /**
@@ -288,9 +447,10 @@ KFR_FUNCTION expression_biquads_l<filters, T, E1> biquad_l(const biquad_params<T
  * @note This implementation has zero latency
  */
 template <size_t filters, typename T, typename E1>
-KFR_FUNCTION expression_biquads<filters, T, E1> biquad(const biquad_params<T> (&bq)[filters], E1&& e1)
+KFR_BIQUAD_DEPRECATED KFR_FUNCTION expression_iir<filters, T, E1> biquad(
+    const biquad_section<T> (&bq)[filters], E1&& e1)
 {
-    return expression_biquads<filters, T, E1>(bq, std::forward<E1>(e1));
+    return expression_iir<filters, T, E1>(std::forward<E1>(e1), iir_state{ iir_params{ bq } });
 }
 
 /**
@@ -300,39 +460,57 @@ KFR_FUNCTION expression_biquads<filters, T, E1> biquad(const biquad_params<T> (&
  * @note This implementation has zero latency
  */
 template <size_t maxfiltercount = 4, typename T, typename E1>
-KFR_FUNCTION expression_handle<T, 1> biquad(const biquad_params<T>* bq, size_t count, E1&& e1)
+KFR_BIQUAD_DEPRECATED KFR_FUNCTION expression_handle<T, 1> biquad(const biquad_section<T>* bq, size_t count,
+                                                                  E1&& e1)
 {
-    constexpr csizes_t<1, 2, 4, 8, 16, 32, 64> sizes;
+    KFR_LOGIC_CHECK(next_poweroftwo(count) <= maxfiltercount,
+                    "biquad: too many biquad sections. Use higher maxfiltercount");
     return cswitch(
-        cfilter(sizes, sizes <= csize_t<maxfiltercount>{}), next_poweroftwo(count),
+        cfilter(internal_generic::biquad_sizes, internal_generic::biquad_sizes <= csize_t<maxfiltercount>{}),
+        next_poweroftwo(count),
         [&](auto x)
         {
             constexpr size_t filters = x;
-            return to_handle(expression_biquads<filters, T, E1>(biquad_block<T, filters>(bq, count),
-                                                                std::forward<E1>(e1)));
+            return to_handle(expression_iir<filters, T, E1>(std::forward<E1>(e1),
+                                                            iir_state{ iir_params<T, filters>(bq, count) }));
         },
         [&] { return to_handle(fixshape(zeros<T>(), fixed_shape<infinite_size>)); });
 }
 
 template <size_t maxfiltercount = 4, typename T, typename E1>
-KFR_FUNCTION expression_handle<T, 1> biquad(const std::vector<biquad_params<T>>& bq, E1&& e1)
+KFR_BIQUAD_DEPRECATED KFR_FUNCTION expression_handle<T, 1> biquad(const std::vector<biquad_section<T>>& bq,
+                                                                  E1&& e1)
 {
     return biquad<maxfiltercount>(bq.data(), bq.size(), std::forward<E1>(e1));
 }
 
+template <size_t filters, typename T, typename E1>
+using expression_biquads_l = expression_iir_l<filters, T, E1>;
+
+template <size_t filters, typename T, typename E1>
+using expression_biquads = expression_iir<filters, T, E1>;
+
 } // namespace CMT_ARCH_NAME
 
 template <typename T>
-class biquad_filter : public expression_filter<T>
+using biquad_params [[deprecated("biquad_params is deprecated. Use biquad_section")]] = biquad_section<T>;
+
+template <typename T, size_t filters = tag_dynamic_vector>
+using biquad_blocks [[deprecated("biquad_blocks is deprecated. Use iir_params")]] = iir_params<T, filters>;
+
+template <typename T>
+class iir_filter : public expression_filter<T>
 {
 public:
-    biquad_filter(const biquad_params<T>* bq, size_t count);
+    iir_filter(const iir_params<T>& params);
 
-    template <size_t N>
-    biquad_filter(const biquad_params<T> (&bq)[N]) : biquad_filter(bq, N)
+    [[deprecated("iir_filter(bq, count) is deprecated. Use iir_filter(iir_params{bq, count})")]] iir_filter(
+        const biquad_section<T>* bq, size_t count)
+        : iir_filter(iir_params<T>(bq, count))
     {
     }
-
-    biquad_filter(const std::vector<biquad_params<T>>& bq) : biquad_filter(bq.data(), bq.size()) {}
 };
+
+template <typename T>
+using biquad_filter [[deprecated("biquad_filter is deprecated. Use iir_filter")]] = iir_filter<T>;
 } // namespace kfr
diff --git a/include/kfr/dsp/dcremove.hpp b/include/kfr/dsp/dcremove.hpp
index e02ad6fd..d017e3a8 100644
--- a/include/kfr/dsp/dcremove.hpp
+++ b/include/kfr/dsp/dcremove.hpp
@@ -32,10 +32,10 @@ namespace kfr
 {
 
 template <typename E1, typename T = flt_type<expression_value_type<E1>>>
-KFR_INTRINSIC expression_biquads<1, T, E1> dcremove(E1&& e1, double cutoff = 0.00025)
+KFR_INTRINSIC expression_iir<1, T, E1> dcremove(E1&& e1, double cutoff = 0.00025)
 {
     const biquad_params<T> bqs[1] = { biquad_highpass(cutoff, 0.5) };
-    return expression_biquads<1, T, E1>(bqs, std::forward<E1>(e1));
+    return expression_iir<1, T, E1>(bqs, std::forward<E1>(e1));
 }
 
 } // namespace kfr
diff --git a/include/kfr/dsp/delay.hpp b/include/kfr/dsp/delay.hpp
index 3253eb15..830be9d3 100644
--- a/include/kfr/dsp/delay.hpp
+++ b/include/kfr/dsp/delay.hpp
@@ -49,14 +49,14 @@ struct delay_state
     {
     }
 
-    mutable univector<T, Tag> data;
-    mutable size_t cursor;
+    univector<T, Tag> data;
+    size_t cursor;
 };
 
 template <typename T>
 struct delay_state<T, 1, 1>
 {
-    mutable T data = T(0);
+    T data = T(0);
 };
 
 template <size_t delay, typename E, bool stateless, univector_tag STag>
@@ -76,8 +76,8 @@ struct expression_delay : expression_with_arguments<E>, public expression_traits
     using T = value_type;
     using expression_with_arguments<E>::expression_with_arguments;
 
-    expression_delay(E&& e, const delay_state<T, delay, STag>& state)
-        : expression_with_arguments<E>(std::forward<E>(e)), state(state)
+    expression_delay(E&& e, state_holder<delay_state<T, delay, STag>, stateless> state)
+        : expression_with_arguments<E>(std::forward<E>(e)), state(std::move(state))
     {
     }
 
@@ -112,7 +112,7 @@ struct expression_delay : expression_with_arguments<E>, public expression_traits
         return concat_and_slice<0, N>(out, in);
     }
 
-    state_holder<const delay_state<T, delay, STag>, stateless> state;
+    mutable state_holder<delay_state<T, delay, STag>, stateless> state;
 };
 
 template <typename E, bool stateless, univector_tag STag>
@@ -146,7 +146,7 @@ struct expression_delay<1, E, stateless, STag> : expression_with_arguments<E>, e
         self.state->data    = in[N - 1];
         return out;
     }
-    state_holder<const delay_state<T, 1, STag>, stateless> state;
+    mutable state_holder<delay_state<T, 1, STag>, stateless> state;
 };
 
 /**
@@ -167,7 +167,7 @@ KFR_INTRINSIC expression_delay<samples, E1, false, samples> delay(E1&& e1)
 
 /**
  * @brief Returns template expression that applies delay to the input (uses ring buffer in state)
- * @param state delay filter state
+ * @param state delay filter state (taken by reference)
  * @param e1 an input expression
  * @code
  * univector<double, 10> v = counter();
@@ -176,12 +176,32 @@ KFR_INTRINSIC expression_delay<samples, E1, false, samples> delay(E1&& e1)
  * @endcode
  */
 template <size_t samples, typename T, typename E1, univector_tag STag>
-KFR_INTRINSIC expression_delay<samples, E1, true, STag> delay(delay_state<T, samples, STag>& state, E1&& e1)
+KFR_INTRINSIC expression_delay<samples, E1, true, STag> delay(
+    E1&& e1, std::reference_wrapper<delay_state<T, samples, STag>> state)
 {
     static_assert(STag == tag_dynamic_vector || (samples >= 1 && samples < 1024), "");
     return expression_delay<samples, E1, true, STag>(std::forward<E1>(e1), state);
 }
 
+/**
+ * @brief Returns template expression that applies delay to the input (uses ring buffer in state)
+ * @param state delay filter state
+ * @param e1 an input expression
+ * @code
+ * univector<double, 10> v = counter();
+ * delay_state<double, 4> state;
+ * auto d = delay(state, v);
+ * @endcode
+ */
+template <size_t samples, typename T, typename E1, univector_tag STag>
+[[deprecated("delay(state, expr) is deprecated. Use delay(expr, std::ref(state))")]] KFR_INTRINSIC
+    expression_delay<samples, E1, true, STag>
+    delay(delay_state<T, samples, STag>& state, E1&& e1)
+{
+    static_assert(STag == tag_dynamic_vector || (samples >= 1 && samples < 1024), "");
+    return expression_delay<samples, E1, true, STag>(std::forward<E1>(e1), std::ref(state));
+}
+
 /**
  * @brief Returns template expression that applies a fractional delay to the input
  * @param e1 an input expression
@@ -193,7 +213,8 @@ KFR_INTRINSIC expression_short_fir<2, T, expression_value_type<E1>, E1> fracdela
     if (CMT_UNLIKELY(delay < 0))
         delay = 0;
     univector<T, 2> taps({ 1 - delay, delay });
-    return expression_short_fir<2, T, expression_value_type<E1>, E1>(std::forward<E1>(e1), taps);
+    return expression_short_fir<2, T, expression_value_type<E1>, E1>(
+        std::forward<E1>(e1), short_fir_state<2, T, expression_value_type<E1>>{ taps });
 }
 } // namespace CMT_ARCH_NAME
 } // namespace kfr
diff --git a/include/kfr/dsp/fir.hpp b/include/kfr/dsp/fir.hpp
index 59f5ab61..d1a8918b 100644
--- a/include/kfr/dsp/fir.hpp
+++ b/include/kfr/dsp/fir.hpp
@@ -57,35 +57,77 @@ struct short_fir_state
     {
     }
     vec<T, tapcount> taps;
-    mutable vec<U, tapcount - 1> delayline;
+    vec<U, tapcount - 1> delayline;
 };
 
+template <typename T>
+struct fir_params
+{
+    univector<T> taps;
+
+    fir_params(const fir_params&)            = default;
+    fir_params(fir_params&&)                 = default;
+    fir_params& operator=(const fir_params&) = default;
+    fir_params& operator=(fir_params&&)      = default;
+
+    fir_params(const T* data, size_t size) : taps(reverse(make_univector(data, size))) {}
+
+    fir_params(univector<T>&& taps) : taps(std::move(taps))
+    {
+        std::reverse(this->taps.begin(), this->taps.end());
+    }
+
+    template <typename Cont, CMT_HAS_DATA_SIZE(Cont)>
+    fir_params(Cont&& taps) : fir_params(std::data(taps), std::size(taps))
+    {
+    }
+};
+
+template <typename Cont>
+fir_params(Cont&&) -> fir_params<std::remove_const_t<container_value_type<Cont>>>;
+
 template <typename T, typename U = T>
 struct fir_state
 {
-    fir_state(const array_ref<const T>& taps)
-        : taps(taps.size()), delayline(taps.size(), U(0)), delayline_cursor(0)
+    fir_state(const fir_state&)            = default;
+    fir_state(fir_state&&)                 = default;
+    fir_state& operator=(const fir_state&) = default;
+    fir_state& operator=(fir_state&&)      = default;
+
+    fir_state(fir_params<T> params)
+        : params(std::move(params)), delayline(this->params.taps.size(), U(0)), delayline_cursor(0)
     {
-        this->taps = reverse(make_univector(taps.data(), taps.size()));
     }
-    univector<T> taps;
-    mutable univector<U> delayline;
-    mutable size_t delayline_cursor;
+    template <typename Cont, CMT_HAS_DATA_SIZE(Cont)>
+    fir_state(Cont&& taps) : params(std::move(taps)), delayline(params.taps.size(), U(0)), delayline_cursor(0)
+    {
+    }
+    template <typename Cont, CMT_HAS_DATA_SIZE(Cont)>
+    void push_delayline(Cont&& state)
+    {
+        delayline.ringbuf_write(delayline_cursor, std::data(state), std::size(state));
+    }
+    fir_params<T> params;
+    univector<U> delayline;
+    size_t delayline_cursor;
 };
 
+template <typename Cont>
+fir_state(Cont&&) -> fir_state<container_value_type<Cont>>;
+
 template <typename U, univector_tag Tag = tag_dynamic_vector>
 struct moving_sum_state
 {
     moving_sum_state() : delayline({ 0 }), head_cursor(0), tail_cursor(1) {}
-    mutable univector<U, Tag> delayline;
-    mutable size_t head_cursor, tail_cursor;
+    univector<U, Tag> delayline;
+    size_t head_cursor, tail_cursor;
 };
 template <typename U>
 struct moving_sum_state<U, tag_dynamic_vector>
 {
     moving_sum_state(size_t sum_length) : delayline(sum_length, U(0)), head_cursor(0), tail_cursor(1) {}
-    mutable univector<U> delayline;
-    mutable size_t head_cursor, tail_cursor;
+    univector<U> delayline;
+    size_t head_cursor, tail_cursor;
 };
 
 inline namespace CMT_ARCH_NAME
@@ -99,8 +141,8 @@ struct expression_short_fir : expression_with_traits<E1>
     static_assert(expression_traits<E1>::dims == 1, "expression_short_fir requires input with dims == 1");
     constexpr static inline bool random_access = false;
 
-    expression_short_fir(E1&& e1, const short_fir_state<tapcount, T, U>& state)
-        : expression_with_traits<E1>(std::forward<E1>(e1)), state(state)
+    expression_short_fir(E1&& e1, state_holder<short_fir_state<tapcount, T, U>, stateless> state)
+        : expression_with_traits<E1>(std::forward<E1>(e1)), state(std::move(state))
     {
     }
 
@@ -120,7 +162,7 @@ struct expression_short_fir : expression_with_traits<E1>
 
         return out;
     }
-    state_holder<const short_fir_state<tapcount, T, U>, stateless> state;
+    mutable state_holder<short_fir_state<tapcount, T, U>, stateless> state;
 };
 
 template <typename T, typename U, typename E1, bool stateless = false>
@@ -131,8 +173,8 @@ struct expression_fir : expression_with_traits<E1>
     static_assert(expression_traits<E1>::dims == 1, "expression_fir requires input with dims == 1");
     constexpr static inline bool random_access = false;
 
-    expression_fir(E1&& e1, const fir_state<T, U>& state)
-        : expression_with_traits<E1>(std::forward<E1>(e1)), state(state)
+    expression_fir(E1&& e1, state_holder<fir_state<T, U>, stateless> state)
+        : expression_with_traits<E1>(std::forward<E1>(e1)), state(std::move(state))
     {
     }
 
@@ -140,7 +182,7 @@ struct expression_fir : expression_with_traits<E1>
     KFR_INTRINSIC friend vec<U, N> get_elements(const expression_fir& self, shape<1> index,
                                                 axis_params<0, N> sh)
     {
-        const size_t tapcount = self.state->taps.size();
+        const size_t tapcount = self.state->params.taps.size();
         const vec<U, N> input = get_elements(self.first(), index, sh);
 
         vec<U, N> output;
@@ -149,17 +191,17 @@ struct expression_fir : expression_with_traits<E1>
         for (size_t i = 0; i < N; i++)
         {
             self.state->delayline.ringbuf_write(cursor, input[i]);
-            U v =
-                dotproduct(self.state->taps.slice(0, tapcount - cursor), self.state->delayline.slice(cursor));
+            U v = dotproduct(self.state->params.taps.slice(0, tapcount - cursor),
+                             self.state->delayline.slice(cursor));
             if (cursor > 0)
-                v = v + dotproduct(self.state->taps.slice(tapcount - cursor),
+                v = v + dotproduct(self.state->params.taps.slice(tapcount - cursor),
                                    self.state->delayline.slice(0, cursor));
             output[i] = v;
         }
         self.state->delayline_cursor = cursor;
         return output;
     }
-    state_holder<const fir_state<T, U>, stateless> state;
+    mutable state_holder<fir_state<T, U>, stateless> state;
 };
 
 template <typename U, typename E1, univector_tag STag, bool stateless = false>
@@ -170,8 +212,8 @@ struct expression_moving_sum : expression_with_traits<E1>
     static_assert(expression_traits<E1>::dims == 1, "expression_moving_sum requires input with dims == 1");
     constexpr static inline bool random_access = false;
 
-    expression_moving_sum(E1&& e1, const moving_sum_state<U, STag>& state)
-        : expression_with_traits<E1>(std::forward<E1>(e1)), state(state)
+    expression_moving_sum(E1&& e1, state_holder<moving_sum_state<U, STag>, stateless> state)
+        : expression_with_traits<E1>(std::forward<E1>(e1)), state(std::move(state))
     {
     }
 
@@ -205,40 +247,94 @@ struct expression_moving_sum : expression_with_traits<E1>
         self.state->tail_cursor = rcursor;
         return output;
     }
-    state_holder<const moving_sum_state<U, STag>, stateless> state;
+    mutable state_holder<moving_sum_state<U, STag>, stateless> state;
 };
 
 /**
  * @brief Returns template expression that applies FIR filter to the input
  * @param e1 an input expression
- * @param taps coefficients for the FIR filter
+ * @param taps coefficients for the FIR filter (taken by value)
  */
-template <typename T, typename E1, univector_tag Tag>
-KFR_INTRINSIC expression_fir<T, expression_value_type<E1>, E1> fir(E1&& e1, const univector<T, Tag>& taps)
+template <typename E1, typename U = expression_value_type<E1>, typename Taps,
+          typename T = std::remove_cv_t<container_value_type<Taps>>>
+[[deprecated("fir(expr, taps) is deprecated. Use fir(expr, fir_params{taps})")]] KFR_INTRINSIC expression_fir<
+    T, U, E1, false>
+fir(E1&& e1, Taps&& taps)
 {
-    return expression_fir<T, expression_value_type<E1>, E1>(std::forward<E1>(e1), taps.ref());
+    return expression_fir<T, U, E1, false>(std::forward<E1>(e1), fir_state<T, U>{ std::forward<Taps>(taps) });
 }
 
 /**
  * @brief Returns template expression that applies FIR filter to the input
- * @param state FIR filter state
  * @param e1 an input expression
+ * @param state coefficients for the FIR filter (taken by value)
  */
-template <typename T, typename U, typename E1>
-KFR_INTRINSIC expression_fir<T, U, E1, true> fir(fir_state<T, U>& state, E1&& e1)
+template <typename T, typename E1, typename U = expression_value_type<E1>>
+KFR_INTRINSIC expression_fir<T, U, E1, false> fir(E1&& e1, fir_params<T> state)
 {
+    return expression_fir<T, U, E1, false>(std::forward<E1>(e1), fir_state<T, U>{ std::move(state) });
+}
+
+/**
+ * @brief Returns template expression that applies FIR filter to the input
+ * @param e1 an input expression
+ * @param state coefficients and state of the filter (taken by reference, ensure proper lifetime)
+ */
+template <typename T, typename E1, typename U>
+KFR_INTRINSIC expression_fir<T, U, E1, true> fir(E1&& e1, std::reference_wrapper<fir_state<T, U>> state)
+{
+    static_assert(std::is_same_v<U, expression_value_type<E1>>, "fir: type mismatch");
     return expression_fir<T, U, E1, true>(std::forward<E1>(e1), state);
 }
 
+/**
+ * @brief Returns template expression that applies FIR filter to the input
+ * @param state FIR filter state (state is referenced, ensure proper lifetime)
+ * @param e1 an input expression
+ */
+template <typename T, typename U, typename E1>
+[[deprecated("fir(state, expr) is deprecated. Use fir(expr, std::ref(state))")]] KFR_INTRINSIC expression_fir<
+    T, U, E1, true>
+fir(fir_state<T, U>& state, E1&& e1)
+{
+    return fir(std::forward<E1>(e1), std::reference_wrapper<fir_state<T, U>>(state));
+}
+
+/**
+ * @brief Returns template expression that performs moving sum on the input
+ * @param e1 an input expression
+ */
+template <typename E1>
+KFR_INTRINSIC expression_moving_sum<expression_value_type<E1>, E1, tag_dynamic_vector> moving_sum(
+    E1&& e1, size_t sum_length)
+{
+    return expression_moving_sum<expression_value_type<E1>, E1, tag_dynamic_vector>(
+        std::forward<E1>(e1), moving_sum_state<expression_value_type<E1>, tag_dynamic_vector>{ sum_length });
+}
+
 /**
  * @brief Returns template expression that performs moving sum on the input
  * @param e1 an input expression
  */
 template <size_t sum_length, typename E1>
-KFR_INTRINSIC expression_moving_sum<expression_value_type<E1>, E1, tag_dynamic_vector> moving_sum(E1&& e1)
+[[deprecated("moving_sum<len> is deprecated. Use moving_sum(expr, len) instead")]] KFR_INTRINSIC
+    expression_moving_sum<expression_value_type<E1>, E1, tag_dynamic_vector>
+    moving_sum(E1&& e1)
+{
+    return expression_moving_sum<expression_value_type<E1>, E1, tag_dynamic_vector>(
+        std::forward<E1>(e1), moving_sum_state<expression_value_type<E1>, tag_dynamic_vector>{ sum_length });
+}
+
+/**
+ * @brief Returns template expression that performs moving sum on the input
+ * @param e1 an input expression
+ * @param state State (taken by reference)
+ */
+template <typename E1, typename U, size_t Tag>
+KFR_INTRINSIC expression_moving_sum<U, E1, Tag, true> moving_sum(
+    E1&& e1, std::reference_wrapper<moving_sum_state<U, Tag>> state)
 {
-    return expression_moving_sum<expression_value_type<E1>, E1, tag_dynamic_vector>(std::forward<E1>(e1),
-                                                                                    sum_length);
+    return expression_moving_sum<expression_value_type<E1>, E1, Tag, true>(std::forward<E1>(e1), state);
 }
 
 /**
@@ -247,9 +343,11 @@ KFR_INTRINSIC expression_moving_sum<expression_value_type<E1>, E1, tag_dynamic_v
  * @param e1 an input expression
  */
 template <typename U, typename E1, univector_tag STag>
-KFR_INTRINSIC expression_moving_sum<U, E1, STag, true> moving_sum(moving_sum_state<U, STag>& state, E1&& e1)
+[[deprecated("moving_sum(state, expr) is deprecated. Use moving_sum(expr, std::ref(state)) "
+             "instead")]] KFR_INTRINSIC expression_moving_sum<U, E1, STag, true>
+moving_sum(moving_sum_state<U, STag>& state, E1&& e1)
 {
-    return expression_moving_sum<U, E1, STag, true>(std::forward<E1>(e1), state);
+    return moving_sum(std::forward<E1>(e1), std::ref(state));
 }
 
 /**
@@ -258,13 +356,28 @@ KFR_INTRINSIC expression_moving_sum<U, E1, STag, true> moving_sum(moving_sum_sta
  * @param e1 an input expression
  * @param taps coefficients for the FIR filter
  */
-template <typename T, size_t TapCount, typename E1>
-KFR_INTRINSIC expression_short_fir<next_poweroftwo(TapCount - 1) + 1, T, expression_value_type<E1>, E1>
-short_fir(E1&& e1, const univector<T, TapCount>& taps)
+template <typename T, size_t TapCount, typename E1,
+          size_t InternalTapCount = next_poweroftwo(TapCount - 1) + 1, typename U = expression_value_type<E1>>
+KFR_INTRINSIC expression_short_fir<InternalTapCount, T, U, E1, false> short_fir(
+    E1&& e1, const univector<T, TapCount>& taps)
 {
     static_assert(TapCount >= 2 && TapCount <= 33, "Use short_fir only for small FIR filters");
-    return expression_short_fir<next_poweroftwo(TapCount - 1) + 1, T, expression_value_type<E1>, E1>(
-        std::forward<E1>(e1), taps);
+    return expression_short_fir<InternalTapCount, T, U, E1, false>(
+        std::forward<E1>(e1), short_fir_state<InternalTapCount, T, U>{ taps });
+}
+/**
+ * @brief Returns template expression that applies FIR filter to the input (count of coefficients must be in
+ * range 2..32)
+ * @param e1 an input expression
+ * @param state FIR filter state (state is referenced, ensure proper lifetime)
+ */
+template <typename T, size_t TapCount, typename E1, typename U>
+KFR_INTRINSIC expression_short_fir<TapCount, T, U, E1, true> short_fir(
+    E1&& e1, std::reference_wrapper<short_fir_state<TapCount, T, U>> state)
+{
+    static_assert(std::is_same_v<U, expression_value_type<E1>>, "short_fir: type mismatch");
+    static_assert(TapCount >= 2 && TapCount <= 33, "Use short_fir only for small FIR filters");
+    return expression_short_fir<TapCount, T, U, E1, true>(std::forward<E1>(e1), state);
 }
 
 /**
@@ -273,13 +386,14 @@ short_fir(E1&& e1, const univector<T, TapCount>& taps)
  * @param state FIR filter state
  * @param e1 an input expression
  */
-template <size_t TapCount, typename T, typename U, typename E1>
-KFR_INTRINSIC expression_short_fir<next_poweroftwo(TapCount - 1) + 1, T, expression_value_type<E1>, E1, true>
-short_fir(short_fir_state<next_poweroftwo(TapCount - 1) + 1, T, U>& state, E1&& e1)
+template <size_t TapCount, size_t InternalTapCount = next_poweroftwo(TapCount - 1) + 1, typename T,
+          typename U, typename E1>
+[[deprecated("short_fir(state, expr) is deprecated, use short_fir(expr, std::ref(state))")]] KFR_INTRINSIC
+    expression_short_fir<InternalTapCount, T, expression_value_type<E1>, E1, true>
+    short_fir(short_fir_state<InternalTapCount, T, U>& state, E1&& e1)
 {
-    static_assert(TapCount >= 2 && TapCount <= 33, "Use short_fir only for small FIR filters");
-    return expression_short_fir<next_poweroftwo(TapCount - 1) + 1, T, expression_value_type<E1>, E1, true>(
-        std::forward<E1>(e1), state);
+    static_assert(InternalTapCount == next_poweroftwo(TapCount - 1) + 1, "short_fir: TapCount mismatch");
+    return short_fir(std::forward<E1>(e1), std::ref(state));
 }
 
 } // namespace CMT_ARCH_NAME
@@ -288,9 +402,10 @@ template <typename T, typename U = T>
 class fir_filter : public filter<U>
 {
 public:
-    fir_filter(const univector_ref<const T>& taps) : state(taps) {}
+    fir_filter(fir_state<T, U> state) : state(std::move(state)) {}
 
-    void set_taps(const univector_ref<const T>& taps) { state = fir_state<T, U>(taps); }
+    void set_taps(fir_params<T> params) { state = std::move(params); }
+    void set_params(fir_params<T> params) { state = std::move(params); }
 
     /// Reset internal filter state
     void reset() final
diff --git a/include/kfr/dsp/iir.hpp b/include/kfr/dsp/iir.hpp
new file mode 100644
index 00000000..2aa0eca3
--- /dev/null
+++ b/include/kfr/dsp/iir.hpp
@@ -0,0 +1,28 @@
+/** @addtogroup biquad
+ *  @{
+ */
+/*
+  Copyright (C) 2016-2023 Dan Cazarin (https://www.kfrlib.com)
+  This file is part of KFR
+
+  KFR is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  KFR is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with KFR.
+
+  If GPL is not suitable for your project, you must purchase a commercial license to use KFR.
+  Buying a commercial license is mandatory as soon as you develop commercial activities without
+  disclosing the source code of your own applications.
+  See https://www.kfrlib.com for details.
+ */
+#pragma once
+
+#include "biquad.hpp"
diff --git a/include/kfr/dsp/iir_design.hpp b/include/kfr/dsp/iir_design.hpp
index 0de49902..86e8a399 100644
--- a/include/kfr/dsp/iir_design.hpp
+++ b/include/kfr/dsp/iir_design.hpp
@@ -1079,7 +1079,7 @@ KFR_FUNCTION zpk<T> iir_bandstop(const zpk<T>& filter, identity<T> lowfreq, iden
 }
 
 template <typename T>
-KFR_FUNCTION std::vector<biquad_params<T>> to_sos(const zpk<T>& filter)
+KFR_FUNCTION iir_params<T> to_sos(const zpk<T>& filter)
 {
     if (filter.p.empty() && filter.z.empty())
         return { biquad_params<T>(filter.k, T(0.), T(0.), T(1.), T(0.), 0) };
@@ -1206,13 +1206,26 @@ KFR_FUNCTION std::vector<biquad_params<T>> to_sos(const zpk<T>& filter)
         pairs[si].z2 = z2;
     }
 
-    std::vector<biquad_params<T>> result(n_sections);
+    iir_params<T> result(n_sections);
     for (size_t si = 0; si < n_sections; si++)
     {
         result[si] = internal::zpk2tf(pairs[n_sections - 1 - si], si == 0 ? filt.k : T(1));
     }
     return result;
 }
+
+/**
+ * @brief Returns template expressions that applies biquad filter to the input.
+ * @param e1 Input expression
+ * @param params IIR filter in ZPK form
+ * @remark This overload converts ZPK to biquad coefficients using to_sos function at every call
+ */
+template <typename T, typename E1>
+KFR_FUNCTION expression_handle<T, 1> iir(E1&& e1, const zpk<T>& params)
+{
+    return iir(std::forward<E1>(e1), to_sos(params));
+}
+
 } // namespace CMT_ARCH_NAME
 
 } // namespace kfr
diff --git a/sources.cmake b/sources.cmake
index 67709216..34ffa1aa 100644
--- a/sources.cmake
+++ b/sources.cmake
@@ -66,6 +66,7 @@ set(
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/fir.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/fir_design.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/goertzel.hpp
+    ${PROJECT_SOURCE_DIR}/include/kfr/dsp/iir.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/iir_design.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/mixdown.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/oscillators.hpp
@@ -262,6 +263,7 @@ set(
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/fir.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/fir_design.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/goertzel.hpp
+    ${PROJECT_SOURCE_DIR}/include/kfr/dsp/iir.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/iir_design.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/mixdown.hpp
     ${PROJECT_SOURCE_DIR}/include/kfr/dsp/oscillators.hpp
diff --git a/src/dsp/biquad.cpp b/src/dsp/biquad.cpp
index 3b310bad..e8f7ed2f 100644
--- a/src/dsp/biquad.cpp
+++ b/src/dsp/biquad.cpp
@@ -31,7 +31,7 @@ namespace kfr
 
 CMT_MULTI_PROTO(namespace impl {
     template <typename T>
-    expression_handle<T, 1> create_biquad_filter(const biquad_params<T>* bq, size_t count);
+    expression_handle<T, 1> create_iir_filter(const iir_params<T>& params);
 } // namespace impl
 )
 
@@ -40,26 +40,25 @@ inline namespace CMT_ARCH_NAME
 namespace impl
 {
 template <typename T>
-expression_handle<T, 1> create_biquad_filter(const biquad_params<T>* bq, size_t count)
+expression_handle<T, 1> create_iir_filter(const iir_params<T>& params)
 {
-    KFR_LOGIC_CHECK(count <= 64, "Too many biquad filters: ", count);
-    return biquad<64>(bq, count, placeholder<T>());
+    return iir(placeholder<T>(), params);
 }
-template expression_handle<float, 1> create_biquad_filter<float>(const biquad_params<float>*, size_t);
-template expression_handle<double, 1> create_biquad_filter<double>(const biquad_params<double>*, size_t);
+template expression_handle<float, 1> create_iir_filter<float>(const iir_params<float>& params);
+template expression_handle<double, 1> create_iir_filter<double>(const iir_params<double>& params);
 } // namespace impl
 } // namespace CMT_ARCH_NAME
 
 #ifdef CMT_MULTI_NEEDS_GATE
 
 template <typename T>
-biquad_filter<T>::biquad_filter(const biquad_params<T>* bq, size_t count)
+iir_filter<T>::iir_filter(const iir_params<T>& params)
 {
-    CMT_MULTI_GATE(this->filter_expr = ns::impl::create_biquad_filter<T>(bq, count));
+    CMT_MULTI_GATE(this->filter_expr = ns::impl::create_iir_filter<T>(params));
 }
 
-template biquad_filter<float>::biquad_filter(const biquad_params<float>*, size_t);
-template biquad_filter<double>::biquad_filter(const biquad_params<double>*, size_t);
+template iir_filter<float>::iir_filter(const iir_params<float>&);
+template iir_filter<double>::iir_filter(const iir_params<double>&);
 
 #endif
 
diff --git a/src/dsp/fir.cpp b/src/dsp/fir.cpp
index 1be6a8d8..0f1d11c1 100644
--- a/src/dsp/fir.cpp
+++ b/src/dsp/fir.cpp
@@ -50,12 +50,12 @@ namespace impl
 template <typename T, typename U>
 void fir_filter<T, U>::process_buffer_impl(U* dest, const U* src, size_t size)
 {
-    make_univector(dest, size) = fir(this->state, make_univector(src, size));
+    make_univector(dest, size) = fir(make_univector(src, size), std::ref(this->state));
 }
 template <typename T, typename U>
 void fir_filter<T, U>::process_expression_impl(U* dest, const expression_handle<U, 1>& src, size_t size)
 {
-    make_univector(dest, size) = fir(this->state, src);
+    make_univector(dest, size) = fir(src, std::ref(this->state));
 }
 
 template class fir_filter<float, float>;
@@ -73,12 +73,12 @@ template class fir_filter<double, complex<double>>;
 template <typename T, typename U>
 void fir_filter<T, U>::process_buffer(U* dest, const U* src, size_t size)
 {
-    make_univector(dest, size) = fir(this->state, make_univector(src, size));
+    make_univector(dest, size) = fir(make_univector(src, size), std::ref(this->state));
 }
 template <typename T, typename U>
 void fir_filter<T, U>::process_expression(U* dest, const expression_handle<U, 1>& src, size_t size)
 {
-    make_univector(dest, size) = fir(this->state, src);
+    make_univector(dest, size) = fir(src, std::ref(this->state));
 }
 template class fir_filter<float, float>;
 template class fir_filter<double, double>;
diff --git a/tests/unit/dsp/biquad.cpp b/tests/unit/dsp/biquad.cpp
index 974c6b59..389c94fe 100644
--- a/tests/unit/dsp/biquad.cpp
+++ b/tests/unit/dsp/biquad.cpp
@@ -66,9 +66,11 @@ TEST(biquad_lowpass1)
                           +0xb.5f265b1be1728p-23, +0xd.d2cb83f8483f8p-24,
                       };
 
-                      const univector<T, size> ir = biquad(bq, unitimpulse<T>());
+                      const univector<T, size> ir  = biquad(bq, unitimpulse<T>());
+                      const univector<T, size> ir2 = iir(unitimpulse<T>(), iir_params{ bq });
 
                       CHECK(absmaxof(choose_array<T>(test_vector_f32, test_vector_f64) - ir) == 0);
+                      CHECK(absmaxof(choose_array<T>(test_vector_f32, test_vector_f64) - ir2) == 0);
                   });
 }
 
diff --git a/tests/unit/dsp/fir.cpp b/tests/unit/dsp/fir.cpp
index f845c847..01237259 100644
--- a/tests/unit/dsp/fir.cpp
+++ b/tests/unit/dsp/fir.cpp
@@ -12,6 +12,27 @@ namespace kfr
 inline namespace CMT_ARCH_NAME
 {
 
+TEST(fir_state)
+{
+    fir_state<float, float> state(univector<float>{ 1, 2, 3, 4, 5, 6 });
+    {
+        const expression_fir<float, float, expression_dimensions<1, float>, false> stateful(
+            dimensions<1>(0.f), state);
+        CHECK(&stateful.state->delayline_cursor != &state.delayline_cursor);
+
+        const expression_fir<float, float, expression_dimensions<1, float>, true> stateless(
+            dimensions<1>(0.f), std::ref(state));
+        CHECK(&stateless.state->delayline_cursor == &state.delayline_cursor);
+    }
+    {
+        auto stateful = fir(dimensions<1>(0.f), state.params);
+        CHECK(&stateful.state->delayline_cursor != &state.delayline_cursor);
+
+        auto stateless = fir(dimensions<1>(0.f), std::ref(state));
+        CHECK(&stateless.state->delayline_cursor == &state.delayline_cursor);
+    }
+}
+
 TEST(fir)
 {
 #ifdef CMT_COMPILER_IS_MSVC