Annotate Artifact Blame Check-ins Using Doc Download Hex

File src/remiaudio/interpolation.cr from the latest check-in

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
   100
   101
   102
   103
   104
   105
   106
   107
   108
   109
   110
   111
   112
   113
   114
   115
   116
   117
   118
   119
   120
   121
   122
   123
   124
   125
   126
   127
   128
   129
   130

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