diff --git a/README.md b/README.md
index df100c0d..dbf4df79 100644
--- a/README.md
+++ b/README.md
@@ -48,6 +48,7 @@ corresponding to `(u,t)` pairs.
 
   - `ConstantInterpolation(u,t)` - A piecewise constant interpolation.
 
+  - `SmoothedConstantInterpolation(u,t)` - An integral preserving continuously differentiable approximation of constant interpolation.
   - `LinearInterpolation(u,t)` - A linear interpolation.
   - `QuadraticInterpolation(u,t)` - A quadratic interpolation.
   - `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`.
diff --git a/docs/src/index.md b/docs/src/index.md
index 1b19d50f..5ce66a6d 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -18,6 +18,7 @@ corresponding to `(u,t)` pairs.
 
   - `ConstantInterpolation(u,t)` - A piecewise constant interpolation.
 
+  - `SmoothedConstantInterpolation(u,t)` - An integral preserving continuously differentiable approximation of constant interpolation.
   - `LinearInterpolation(u,t)` - A linear interpolation.
   - `QuadraticInterpolation(u,t)` - A quadratic interpolation.
   - `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`.
diff --git a/docs/src/manual.md b/docs/src/manual.md
index 3818a9e3..e599fad9 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -6,6 +6,7 @@ QuadraticInterpolation
 LagrangeInterpolation
 AkimaInterpolation
 ConstantInterpolation
+SmoothedConstantInterpolation
 QuadraticSpline
 CubicSpline
 BSplineInterpolation
diff --git a/docs/src/methods.md b/docs/src/methods.md
index 648131eb..9a7b54ff 100644
--- a/docs/src/methods.md
+++ b/docs/src/methods.md
@@ -77,6 +77,24 @@ A = ConstantInterpolation(u, t, dir = :right)
 plot(A)
 ```
 
+## Smoothed Constant Interpolation
+
+This function is much like the constant interpolation above, but the transition
+between consecutive values is smoothed out so that the function is continuously
+differentiable. The smoothing is done in such a way that the integral of this function
+is never much off from the same integral of constant interpolation without smoothing (because of the symmetry of the smoothing sections).
+The maximum smoothing distance in the `t` direction from the data points can be set
+with the keyword argument `d_max`.
+
+```@example tutorial
+A = ConstantInterpolation(u, t)
+plot(A)
+A = SmoothedConstantInterpolation(u, t; d_max = 10)
+plot!(A)
+```
+
+Note that `u[end]` is ignored.
+
 ## Quadratic Spline
 
 This is the quadratic spline. It is a continuously differentiable interpolation
diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl
index 22dc7345..d25e8a17 100644
--- a/src/DataInterpolations.jl
+++ b/src/DataInterpolations.jl
@@ -101,10 +101,10 @@ function Base.showerror(io::IO, ::IntegralNotInvertibleError)
 end
 
 export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
-       AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
-       BSplineInterpolation, BSplineApprox, CubicHermiteSpline, PCHIPInterpolation,
-       QuinticHermiteSpline, LinearInterpolationIntInv, ConstantInterpolationIntInv,
-       ExtrapolationType
+       AkimaInterpolation, ConstantInterpolation, SmoothedConstantInterpolation,
+       QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox,
+       CubicHermiteSpline, PCHIPInterpolation, QuinticHermiteSpline,
+       LinearInterpolationIntInv, ConstantInterpolationIntInv, ExtrapolationType
 
 # added for RegularizationSmooth, JJS 11/27/21
 ### Regularization data smoothing and interpolation
diff --git a/src/derivatives.jl b/src/derivatives.jl
index a81ef368..3439aa27 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -67,6 +67,19 @@ function _derivative(A::LinearInterpolation, t::Number, iguess)
     slope
 end
 
+function _derivative(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Number, iguess)
+    idx = get_idx(A, t, iguess)
+    d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx)
+
+    if (t - A.t[idx]) < d_lower
+        -2c_lower * ((t - A.t[idx]) / d_lower - 1) / d_lower
+    elseif (A.t[idx + 1] - t) < d_upper
+        2c_upper * (1 - (A.t[idx + 1] - t) / d_upper) / d_upper
+    else
+        zero(c_upper / oneunit(t))
+    end
+end
+
 function _derivative(A::QuadraticInterpolation, t::Number, iguess)
     idx = get_idx(A, t, iguess)
     Δt = t - A.t[idx]
diff --git a/src/integrals.jl b/src/integrals.jl
index 41e52e4a..6019867f 100644
--- a/src/integrals.jl
+++ b/src/integrals.jl
@@ -152,6 +152,38 @@ function _integral(
     end
 end
 
+function _integral(A::SmoothedConstantInterpolation{<:AbstractVector},
+        idx::Number, t1::Number, t2::Number)
+    d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx)
+
+    bound_lower = A.t[idx] + d_lower
+    bound_upper = A.t[idx + 1] - d_upper
+
+    out = A.u[idx] * (t2 - t1)
+
+    # Fix extrapolation behavior as constant for now
+    if t1 <= first(A.t)
+        t1 = first(A.t)
+    elseif t2 >= last(A.t)
+        t2 = last(A.t)
+    end
+
+    if t1 < bound_lower
+        t2_ = min(t2, bound_lower)
+        out -= c_lower * d_lower *
+               (((t2_ - A.t[idx]) / d_lower - 1)^3 - ((t1 - A.t[idx]) / d_lower - 1)^3) / 3
+    end
+
+    if t2 > bound_upper
+        t1_ = max(t1, bound_upper)
+        out += c_upper * d_upper *
+               ((1 - (A.t[idx + 1] - t2) / d_upper)^3 -
+                (1 - (A.t[idx + 1] - t1_) / d_upper)^3) / 3
+    end
+
+    out
+end
+
 function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}},
         idx::Number, t1::Number, t2::Number)
     α, β = get_parameters(A, idx)
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index ad2e1eb0..28a44bab 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -366,6 +366,79 @@ function ConstantInterpolation(
         cache_parameters, assume_linear_t)
 end
 
+"""
+    SmoothedConstantInterpolation(u, t; d_max = Inf, extrapolate = false,
+        cache_parameters = false, assume_linear_t = 1e-2)
+
+It is a method for interpolating constantly with forward fill, with smoothing around the
+value transitions to make the curve continuously differentiable while the integral never
+drifts far from the integral of constant interpolation. In this interpolation type, `u[end]` is ignored.
+
+## Arguments
+
+  - `u`: data points.
+  - `t`: time points.
+
+## Keyword Arguments
+
+  - `d_max`: Around each time point `tᵢ` there is a continuously differentiable (quadratic) transition between `uᵢ₋₁` and `uᵢ`,
+    on the interval `[tᵢ - d, tᵢ + d]`. The distance `d` is determined as `d = min((tᵢ - tᵢ₋₁)/2, (tᵢ₊₁ - tᵢ)/2, d_max)`.
+  - `extrapolation`: The extrapolation type applied left and right of the data. Possible options
+    are `ExtrapolationType.None` (default), `ExtrapolationType.Constant`, `ExtrapolationType.Linear`
+    `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` (also made smooth at the boundaries) and `ExtrapolationType.Reflective`.
+  - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for
+    the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`.
+  - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for
+    the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`.
+  - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behavior for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
+"""
+struct SmoothedConstantInterpolation{uType, tType, IType, dType, cType, dmaxType, T, N} <:
+       AbstractInterpolation{T, N}
+    u::uType
+    t::tType
+    I::IType
+    p::SmoothedConstantParameterCache{dType, cType}
+    d_max::dmaxType
+    extrapolation_left::ExtrapolationType.T
+    extrapolation_right::ExtrapolationType.T
+    iguesser::Guesser{tType}
+    cache_parameters::Bool
+    linear_lookup::Bool
+    function SmoothedConstantInterpolation(
+            u, t, I, p, d_max, extrapolation_left,
+            extrapolation_right, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(I), typeof(p.d),
+            typeof(p.c), typeof(d_max), eltype(u), N}(
+            u, t, I, p, d_max, extrapolation_left, extrapolation_right,
+            Guesser(t), cache_parameters, linear_lookup)
+    end
+end
+
+function SmoothedConstantInterpolation(
+        u, t; d_max = Inf, extrapolation::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
+        cache_parameters = false, assume_linear_t = 1e-2)
+    extrapolation_left, extrapolation_right = munge_extrapolation(
+        extrapolation, extrapolation_left, extrapolation_right)
+    u, t = munge_data(u, t)
+    p = SmoothedConstantParameterCache(
+        u, t, cache_parameters, d_max, extrapolation_left, extrapolation_right)
+    A = SmoothedConstantInterpolation(
+        u, t, nothing, p, d_max, extrapolation_left,
+        extrapolation_right, cache_parameters, assume_linear_t)
+    I = cumulative_integral(A, cache_parameters)
+    SmoothedConstantInterpolation(
+        u, t, I, p, d_max, extrapolation_left,
+        extrapolation_right, cache_parameters, assume_linear_t)
+end
+
 """
     QuadraticSpline(u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false)
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index 14467311..ccbc668f 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -164,7 +164,7 @@ function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess
     @evalpoly wj A.u[idx] A.b[idx] A.c[idx] A.d[idx]
 end
 
-# ConstantInterpolation Interpolation
+# Constant Interpolation
 function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
@@ -187,6 +187,22 @@ function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, igu
     A.u[:, idx]
 end
 
+# Smoothed constant Interpolation
+function _interpolate(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Number, iguess)
+    idx = get_idx(A, t, iguess)
+    d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx)
+
+    out = A.u[idx]
+
+    if (t - A.t[idx]) < d_lower
+        out -= c_lower * ((t - A.t[idx]) / d_lower - 1)^2
+    elseif (A.t[idx + 1] - t) < d_upper
+        out += c_upper * (1 - (A.t[idx + 1] - t) / d_upper)^2
+    end
+
+    out
+end
+
 # QuadraticSpline Interpolation
 function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
     idx = get_idx(A, t, iguess)
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 1a0aaef4..eff36bb8 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -207,6 +207,22 @@ function get_parameters(A::LinearInterpolation, idx)
     end
 end
 
+function get_parameters(A::SmoothedConstantInterpolation, idx)
+    if A.cache_parameters
+        d_lower = A.p.d[idx]
+        d_upper = A.p.d[idx + 1]
+        c_lower = A.p.c[idx]
+        c_upper = A.p.c[idx + 1]
+        d_lower, d_upper, c_lower, c_upper
+    else
+        d_lower, c_lower = smoothed_constant_interpolation_parameters(
+            A.u, A.t, A.d_max, idx, A.extrapolation_left, A.extrapolation_right)
+        d_upper, c_upper = smoothed_constant_interpolation_parameters(
+            A.u, A.t, A.d_max, idx + 1, A.extrapolation_left, A.extrapolation_right)
+        d_lower, d_upper, c_lower, c_upper
+    end
+end
+
 function get_parameters(A::QuadraticInterpolation, idx)
     if A.cache_parameters
         A.p.α[idx], A.p.β[idx]
diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl
index 737f736b..ca1f6c00 100644
--- a/src/parameter_caches.jl
+++ b/src/parameter_caches.jl
@@ -32,6 +32,39 @@ function linear_interpolation_parameters(u::AbstractArray{T, N}, t, idx) where {
     return slope
 end
 
+struct SmoothedConstantParameterCache{dType, cType}
+    d::dType
+    c::cType
+end
+
+function SmoothedConstantParameterCache(
+        u, t, cache_parameters, d_max, extrapolation_left, extrapolation_right)
+    if cache_parameters
+        parameters = smoothed_constant_interpolation_parameters.(
+            Ref(u), Ref(t), d_max, eachindex(t), extrapolation_left, extrapolation_right)
+        d, c = collect.(eachrow(stack(collect.(parameters))))
+        SmoothedConstantParameterCache(d, c)
+    else
+        SmoothedConstantParameterCache(eltype(t)[], eltype(u)[])
+    end
+end
+
+function smoothed_constant_interpolation_parameters(
+        u, t, d_max, idx, extrapolation_left, extrapolation_right)
+    if isone(idx) || (idx == length(t))
+        # If extrapolation is periodic, make the transition differentiable
+        if extrapolation_left == extrapolation_right == ExtrapolationType.Periodic
+            min(t[end] - t[end - 1], t[2] - t[1], 2d_max) / 2, (u[1] - u[end - 1]) / 2
+        else
+            d = isone(idx) ? min(t[2] - t[1], 2d_max) / 2 :
+                min(t[end] - t[end - 1], 2d_max) / 2
+            d, zero(one(eltype(u)) / 2)
+        end
+    else
+        min(t[idx] - t[idx - 1], t[idx + 1] - t[idx], 2d_max) / 2, (u[idx] - u[idx - 1]) / 2
+    end
+end
+
 struct QuadraticParameterCache{pType}
     α::pType
     β::pType
diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl
index 376b3310..3190e601 100644
--- a/test/derivative_tests.jl
+++ b/test/derivative_tests.jl
@@ -31,7 +31,7 @@ function test_derivatives(method; args = [], kwargs = [], name::String)
         # Interpolation time points
         for _t in t[2:(end - 1)]
             if func isa BSplineInterpolation || func isa BSplineApprox ||
-               func isa CubicHermiteSpline
+               func isa CubicHermiteSpline || func isa SmoothedConstantInterpolation
                 fdiff = forward_fdm(5, 1; geom = true)(func, _t)
                 fdiff2 = forward_fdm(5, 1; geom = true)(t -> derivative(func, t), _t)
             else
@@ -150,6 +150,13 @@ end
     @test all(derivative.(Ref(A), t2 .+ 0.1) .== 0.0)
 end
 
+@testset "SmoothedConstantInterpolation" begin
+    u = [5.5, 2.7, 5.1, 3.0]
+    t = [2.5, 5.6, 6.3, 8.9]
+    test_derivatives(SmoothedConstantInterpolation; args = [u, t],
+        name = "Smoothed constant interpolation")
+end
+
 @testset "Quadratic Spline" begin
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
diff --git a/test/integral_tests.jl b/test/integral_tests.jl
index 283b8dce..a528ed66 100644
--- a/test/integral_tests.jl
+++ b/test/integral_tests.jl
@@ -67,6 +67,20 @@ end
         name = "Linear Interpolation (Vector) with random points")
 end
 
+@testset "SmoothedConstantInterpolation" begin
+    u = [1.0, 4.0, 9.0, 16.0]
+    t = [1.0, 2.0, 3.0, 4.0]
+    test_integral(
+        SmoothedConstantInterpolation; args = [u, t], name = "Smoothed constant interpolation")
+
+    A_constant = ConstantInterpolation(u, t)
+    I_ref = DataInterpolations.integral(A_constant, first(t), last(t))
+    I_smoothed = [DataInterpolations.integral(
+                      SmoothedConstantInterpolation(u, t; d_max), first(t), last(t))
+                  for d_max in 0.0:0.1:1.0]
+    @test all(I_smoothed .≈ I_ref)
+end
+
 @testset "QuadraticInterpolation" begin
     u = [1.0, 4.0, 9.0, 16.0]
     t = [1.0, 2.0, 3.0, 4.0]
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index ab0819d5..99afa59e 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -522,6 +522,19 @@ end
     @test A(-Inf) == first(u)
 end
 
+@testset "Smoothed constant Interpolation" begin
+    test_interpolation_type(SmoothedConstantInterpolation)
+    u = [0.0, 2.0, 1.0, 3.0]
+    t = [1.2, 2.5, 5.7, 8.7]
+    d_max = 0.5
+    A = SmoothedConstantInterpolation(u, t; d_max)
+    test_cached_index(A)
+
+    @test A(1.9) == u[1]
+    @test A(3.1) == u[2]
+    @test A(2.5) ≈ (u[1] + u[2]) / 2
+end
+
 @testset "QuadraticSpline Interpolation" begin
     test_interpolation_type(QuadraticSpline)
 
diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl
index 1f8ff574..e7aa6b63 100644
--- a/test/parameter_tests.jl
+++ b/test/parameter_tests.jl
@@ -15,6 +15,14 @@ end
     test_cached_integration(LinearInterpolation, u, t)
 end
 
+@testset "Smoothed constant Interpolation" begin
+    u = [1.0, 5.0, 3.0, 4.0, 4.0]
+    t = collect(1:5)
+    A = SmoothedConstantInterpolation(u, t; cache_parameters = true)
+    @test A.p.d ≈ [0.5, 0.5, 0.5, 0.5, 0.5]
+    @test A.p.c ≈ [0.0, 2.0, -1.0, 0.5, 0.0]
+end
+
 @testset "Quadratic Interpolation" begin
     u = [1.0, 5.0, 3.0, 4.0, 4.0]
     t = collect(1:5)