#### RemiAudio
#### Copyright (C) 2022-2024 Remilia Scarlet <remilia@posteo.jp>
####
#### This program is free software: you can redistribute it and/or modify it
#### under the terms of the GNU Affero General Public License as published by
#### the Free Software Foundation, either version 3 of the License, or (at your
#### option) any later version.
####
#### This program 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 Affero General Public
#### License for more details.
####
#### You should have received a copy of the GNU Affero General Public License
#### along with this program. If not, see <https://www.gnu.org/licenses/>.
module RemiAudio::Interpolate
extend self
# An interpolation type.
enum Mode
Linear
Cubic
Hermite
HermiteAlt
BSpline
Parabolic2X
Optimal2X4P2O
Optimal32X4P4O
end
# Performs linear interpolation.
@[AlwaysInline]
def linear(x1 : Float64, x2 : Float64, mu : Float64) : Float64
(x1 + (mu * (x2 - x1)))
end
# Performs cubic interpolation.
@[AlwaysInline]
def cubic(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
c1 = (x2 - x0) * 0.5
c3 = (x1 - x2) * 1.5 + (x3 - x0) * 0.5
c2 = c1 + x0 - x1 - c3
((c3 * mu + c2) * mu + c1) * mu + x1
end
# Performs hermite interpolation.
#
# This is an alternate version with slightly different coefficients.
@[AlwaysInline]
def hermiteAlt(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
mu2 = mu * mu
c0 = 3 * x1 - 3 * x2 + x3 - x0
c1 = 2 * x0 - 5 * x1 + 4 * x2 - x3
c2 = x2 - x0
c3 = 2 * x1
(c0 * mu * mu2 + c1 * mu2 + c2 * mu + c3) / 2
end
### These below are taken from the "Elephant Paper".
###
### http://yehar.com/blog/wp-content/uploads/2009/08/deip.pdf
# Performs B-Spline interpolation.
@[AlwaysInline]
def bspline(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
ym1py1 = x0 + x2
c0 = 1 / 6.0 * ym1py1 + 2 / 3.0 * x1
c1 = 1 / 2.0 * (x2 - x0)
c2 = 1 / 2.0 * ym1py1 - x1
c3 = 1 / 2.0 * (x1 - x2) + 1 / 6.0 * (x3 - x0)
((c3 * mu + c2) * mu + c1) * mu + c0
end
# Performs hermite interpolation.
@[AlwaysInline]
def hermite(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
# Goddess do I miss prefix notation
c0 = x1
c1 = 0.5 * (x2 - x0)
c2 = x0 - (2.5 * x1) + (2 * x2) - (0.5 * x3)
c3 = (0.5 * (x3 - x0)) + (1.5 * (x1 - x2))
(((((c3 * mu) + c2) * mu) + c1) * mu) + c0
end
# Performs parabolic interpolation. This expects 2x oversampled data.
@[AlwaysInline]
def parabolic2x(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
y1mym1 : Float64 = x2 - x0
c0 : Float64 = (0.5 * x1) + (0.25 * (x0 + x2))
c1 : Float64 = 0.5 * y1mym1
c2 : Float64 = 0.25 * (x1 - x1 - y1mym1)
(((c2 * mu) + c1) * mu) + c0
end
# Performs optimal 4-point 2nd order interpolation. This expects 2x
# oversampled data.
@[AlwaysInline]
def optimal2x4P2O(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
z : Float64 = mu - 0.5
even1 : Float64 = x2 + x1
odd1 : Float64 = x2 - x1
even2 : Float64 = x3 + x0
odd2 : Float64 = x3 - x0
c0 : Float64 = (even1 * 0.42334633257225274_f64) + (even2 * 0.07668732202139628_f64)
c1 : Float64 = (odd1 * 0.26126047291143606_f64) + (odd2 * 0.24778879018226652_f64)
c2 : Float64 = (even1 * -0.213439787561776841_f64) + (even2 * 0.21303593243799016_f64)
(((c2 * z) + c1) * z) + c0
end
# Performs optimal 4-point 4th order interpolation. This expects 32x
# oversampled data.
@[AlwaysInline]
def optimal32x4P4O(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
z : Float64 = mu - 0.5
even1 : Float64 = x2 + x1
odd1 : Float64 = x2 - x1
even2 : Float64 = x3 + x0
odd2 : Float64 = x3 - x0
c0 = (even1 * 0.46835497211269561_f64) + (even2 * 0.03164502784253309_f64)
c1 = (odd1 * 0.56001293337091440_f64) + (odd2 * 0.14666238593949288_f64)
c2 = (even1 * -0.250038759826233691_f64) + (even2 * 0.25003876124297131_f64)
c3 = (odd1 * -0.49949850957839148_f64) + (odd2 * 0.16649935475113800_f64)
c4 = (even1 * 0.00016095224137360_f64) + (even2 * -0.00016095810460478_f64)
(((c4 * z + c3) * z + c2) * z + c1) * z + c0
end
end