RemiAudio: sinc.cr at tip

Annotate Artifact Blame Check-ins Using Doc Download Hex
File src/remiaudio/resampler/sinc.cr from the latest check-in

#### Based on libsamplerate
####
#### Copyright (c) 2002-2021, Erik de Castro Lopo <erikd@mega-nerd.com>
#### Copyright (c) 2024, Remilia Scarlet <remilia@posteo.jp>
#### All rights reserved.
####
#### This code is released under 2-clause BSD license. Please see the
#### file at : https://github.com/libsndfile/libsamplerate/blob/master/COPYING
require "./resampler"

@[Link(ldflags: "#{__DIR__}/libsinctables.a")]
lib LibRASincTables
  $sinc_fast_coeffs_size : LibC::Int
  $sinc_medium_coeffs_size : LibC::Int
  $sinc_best_coeffs_size : LibC::Int

  $sinc_fast_coeffs_index_inc : LibC::Int
  $sinc_medium_coeffs_index_inc : LibC::Int
  $sinc_best_coeffs_index_inc : LibC::Int

  $sinc_fast_coeffs : StaticArray(Float32, 2464)
  $sinc_medium_coeffs : StaticArray(Float32, 22438)
  $sinc_best_coeffs : StaticArray(Float32, 340239)
end

module RemiAudio::Resampler
  # A bandlimited interpolator derived from the mathematical sinc function.
  # Based on the techniques of [Julius
  # O. Smith](http://ccrma.stanford.edu/~jos/resample/), although this code was
  # developed independently.
  abstract class SincResampler < Resampler
    MAX_CHANNELS = 128

    # :nodoc:
    SHIFT_BITS = 12

    # :nodoc:
    FP_ONE = (1 << SHIFT_BITS).to_f64!

    # :nodoc
    INV_FP_ONE = 1.0 / FP_ONE

    # Defines the quality of a `SincResampler`.
    enum Quality
      # The fastest, but lowest quality, for a `SincResampler`. The resampler
      # will have a signal-to-noise ratio of 97dB and a bandwidth of 80%.
      Fast

      # A fairly fast, medium quality, for a `SincResampler`. The resampler will
      # have a signal-to-noise ratio of 97dB and a bandwidth of 90%
      Medium

      # A very slow, but very high quality, for a `SincResampler`.  This
      # provides a worst case signal-to-noise ratio of 97dB at a bandwidth of
      # 97%
      Best
    end

    # :nodoc:
    alias Coefficients = Pointer(Float32)

    @inCount : Int64 = 0i64
    @inUsed : Int64 = 0i64
    @outCount : Int64 = 0i64
    @outGen : Int64 = 0i64
    @coeffHalfLen : Int32 = 0
    @indexInc : Int32 = 0
    @srcRatio : Float64 = 0.0
    @inputIndex : Float64 = 0.0
    @coeffs : Coefficients
    @bCurrent : Int32 = 0
    @bEnd : Int32 = 0
    @bRealEnd : Int32 = 0
    @bLen : Int32 = 0
    @leftCalc : Array(Float64) = Array(Float64).new(MAX_CHANNELS, 0.0)
    @rightCalc : Array(Float64) = Array(Float64).new(MAX_CHANNELS, 0.0)
    @buffer : Array(Float32)

    # The `Quality` for this instance.
    getter quality : Quality

    # Appease the compiler deities
    private def initialize(@channels, @quality, @coeffs, @buffer)
    end

    # :inherit:
    def reset : Nil
      super
      @bCurrent = 0
      @bEnd = 0
      @bRealEnd = -1
      @srcRatio = 0.0
      @inputIndex = 0.0
      @buffer.fill(0.0f32)
    end

    @[AlwaysInline]
    protected def toFp(x : Float64) : Int32
      (x.round * FP_ONE).to_i32!
    end

    @[AlwaysInline]
    protected def toFp(x : Int) : Int32
      (x << SHIFT_BITS).to_i32!
    end

    @[AlwaysInline]
    protected def fpToInt(x : Int32) : Int32
      x >> SHIFT_BITS
    end

    @[AlwaysInline]
    protected def fpFracPart(x : Int32)
      x & ((1 << SHIFT_BITS) &- 1)
    end

    @[AlwaysInline]
    protected def fpToF64(x : Int32) : Float64
      fpFracPart(x) * INV_FP_ONE
    end

    @[AlwaysInline]
    protected def intDivCeil(divident : Int32, divisor : Int32) # == (divident.to_f32! / divisor).ceiling.to_i32!
      RemiLib.assert(divident >= 0 && divisor > 0) # For positive numbers only
      (divident &+ (divisor &- 1)).tdiv(divisor)
    end

    protected def prepareData(data : Data, halfFilterChanLen : Int32) : Nil
      return if @bRealEnd >= 0 # Should be terminating, just return
      return if data.dataIn.empty?

      len : Int32 = 0
      if @bCurrent == 0
        # Initial state. Set up zeros at the start of the buffer and then load
        # new data after that.
        len = @bLen - 2 * halfFilterChanLen
        @bCurrent = halfFilterChanLen
        @bEnd = halfFilterChanLen
      elsif (@bEnd + halfFilterChanLen + @channels) < @bLen
        # Load data at current end position.
        len = Math.max(@bLen - @bCurrent - halfFilterChanLen, 0)
      else
        # Move data at end of buffer back to the start of the buffer.
        len = @bEnd - @bCurrent
        @buffer.to_unsafe.move_from(@buffer.to_unsafe + @bCurrent - halfFilterChanLen,
                                    halfFilterChanLen + len)
        @bCurrent = halfFilterChanLen
        @bEnd = @bCurrent + len

        # Now load data at current end of buffer.
        len = Math.max(@bLen - @bCurrent - halfFilterChanLen, 0)
      end

      len = Math.min((@inCount - @inUsed).to_i32!, len)
      len -= (len % @channels)

      if len < 0 || (@bEnd + len) > @bLen
        raise Error.new("prepareData: bad length")
      end

      (@buffer.to_unsafe + @bEnd).copy_from(data.dataIn.to_unsafe + @inUsed, len)

      @bEnd += len
      @inUsed += len

      if @inUsed == @inCount && (@bEnd - @bCurrent) < (2 * halfFilterChanLen) && data.endOfInput?
        # Handle the case where all data in the current buffer has been consumed
        # and this is the last buffer.

        if (@bLen - @bEnd) < (halfFilterChanLen + 5)
          # If necessary, move data down to the start of the buffer.
          len = @bEnd - @bCurrent
          @buffer.to_unsafe.move_from(@buffer.to_unsafe + @bCurrent - halfFilterChanLen,
                                      halfFilterChanLen + len)

          @bCurrent = halfFilterChanLen
          @bEnd = @bCurrent + len
        end

        @bRealEnd = @bEnd
        len = halfFilterChanLen + 5

        if len < 0 || (@bEnd + len) > @bLen
          len = @bLen - @bEnd
        end

        @buffer.fill(0.0, @bEnd, len)
        @bEnd += len
      end
    end
  end

  # A `SincResampler` optimized for two channels.
  class SincResamplerStereo < SincResampler
    # Creates a new `SincResamplerStereo` instance that will operate at the
    # given `Quality`.
    def initialize(@quality : Quality)
      @channels = 2
      case @quality
      in .fast?
        @coeffs = LibRASincTables.sinc_fast_coeffs.to_unsafe
        @coeffHalfLen = LibRASincTables.sinc_fast_coeffs_size.to_i32 - 2
        @indexInc = LibRASincTables.sinc_fast_coeffs_index_inc.to_i32
      in .medium?
        @coeffs = LibRASincTables.sinc_medium_coeffs.to_unsafe
        @coeffHalfLen = LibRASincTables.sinc_medium_coeffs_size.to_i32 - 2
        @indexInc = LibRASincTables.sinc_medium_coeffs_index_inc.to_i32
      in .best?
        @coeffs = LibRASincTables.sinc_best_coeffs.to_unsafe
        @coeffHalfLen = LibRASincTables.sinc_best_coeffs_size.to_i32 - 2
        @indexInc = LibRASincTables.sinc_best_coeffs_index_inc.to_i32
      end

      @bLen = 3 * ((@coeffHalfLen + 2.0) / @indexInc * SRC_MAX_RATIO + 1).round.to_i32!
      @bLen = Math.max(@bLen, 4096)
      @bLen *= @channels
      @bLen += 1

      @buffer = Array(Float32).new(@bLen * @channels, 0.0f32)

      reset
    end

    private def calcOutput(increment : Int32, startFilterIndex : Int32, scale : Float64, output : Slice(Float32)) : Nil
      fraction : Float64 = 0.0
      left : StaticArray(Float64, 2) = StaticArray(Float64, 2).new(0.0)
      right : StaticArray(Float64, 2) = StaticArray(Float64, 2).new(0.0)
      icoeff : Float64 = 0.0
      idx : Int32 = 0

      # Convert input parameters into fixed point.
      maxFilterIndex : Int32 = toFp(@coeffHalfLen)

      # First apply the left half of the filter.
      filterIndex : Int32 = startFilterIndex
      coeffCount : Int32 = (maxFilterIndex - filterIndex).tdiv(increment)
      filterIndex = filterIndex + coeffCount * increment
      dataIndex : Int32 = @bCurrent - @channels * coeffCount
      if dataIndex < 0 # Avoid underflow access to filter->buffer.
        steps : Int32 = intDivCeil(-dataIndex, 2)
        # If the assert triggers we would have to take care not to underflow/overflow
        RemiLib.assert(steps <= intDivCeil(filterIndex, increment))
        filterIndex -= (increment * steps)
        dataIndex += (steps * 2)
      end

      while filterIndex >= 0
        fraction = fpToF64(filterIndex)
        idx = fpToInt(filterIndex)
        #RemiLib.assert(idx >= 0 && (idx + 1) < (@coeffHalfLen + 2))
        icoeff = @coeffs[idx].to_f64! + fraction * (@coeffs[idx + 1] - @coeffs[idx])
        #RemiLib.assert(dataIndex >= 0 && (dataIndex + 1) < @bLen)
        #RemiLib.assert((dataIndex + 1) < @bEnd)
        2.times do |ch|
          left[ch] += (icoeff * @buffer[dataIndex + ch])
        end

        filterIndex -= increment
        dataIndex = dataIndex + 2
      end

      # Now apply the right half of the filter.
      filterIndex = increment - startFilterIndex
      coeffCount = (maxFilterIndex - filterIndex).tdiv(increment)
      filterIndex = filterIndex + coeffCount * increment
      dataIndex = @bCurrent + @channels * (1 + coeffCount)

      loop do
        fraction = fpToF64(filterIndex)
        idx = fpToInt(filterIndex)
        #RemiLib.assert(idx >= 0 && (idx + 1) < (@coeffHalfLen + 2))
        icoeff = @coeffs[idx].to_f64! + fraction * (@coeffs[idx + 1] - @coeffs[idx])
        #RemiLib.assert(dataIndex >= 0 && (dataIndex + 1) < @bLen)
        #RemiLib.assert("dataIndex = #{dataIndex}, bEnd: #{@bEnd}",
        #               (dataIndex + 1) < @bEnd)
        2.times do |ch|
          right[ch] += (icoeff * @buffer[dataIndex + ch])
        end

        filterIndex -= increment
        dataIndex = dataIndex - 2

        break unless filterIndex > 0
      end

      2.times do |ch|
        output[ch] = (scale * (left[ch] + right[ch])).to_f32!
      end
    end

    protected def variProcess(data : Data) : Nil
      samplesInHand : Int32 = 0
      floatIncrement : Float64 = 0.0
      startFilterIndex : Int32 = 0
      increment : Int32 = 0
      @inCount = data.inputFrames * @channels
      @outCount = data.outputFrames * @channels
      @inUsed = 0
      @outGen = 0

      ratio : Float64 = @lastRatio

      # Check the sample rate ratio wrt the buffer len.
      count : Float64 = (@coeffHalfLen + 2.0) / @indexInc
      if Math.min(@lastRatio, data.srcRatio) < 1.0
        count /= Math.min(@lastRatio, data.srcRatio)
      end

      # Maximum coefficients on either side of center point.
      halfFilterChanLen : Int32 = @channels * (count.round.to_i32! &+ 1)

      inputIndex : Float64 = @lastPosition

      rem : Float64 = fmodOne(inputIndex)
      @bCurrent = ((@bCurrent.to_i64! + @channels * (inputIndex - rem).round.to_i64!) % @bLen).to_i32!
      inputIndex = rem

      terminate : Float64 = 1.0 / ratio + 1e-20

      # Main processing loop.
      while @outGen < @outCount
        # Need to reload buffer?
        samplesInHand = (@bEnd - @bCurrent + @bLen) % @bLen

        if samplesInHand <= halfFilterChanLen
          prepareData(data, halfFilterChanLen)
          samplesInHand = (@bEnd - @bCurrent + @bLen) % @bLen
          break if samplesInHand <= halfFilterChanLen
        end

        # This is the termination condition.
        if @bRealEnd >= 0
          break if (@bCurrent + inputIndex + terminate) >= @bRealEnd
        end

        if @outCount > 0 && (@lastRatio - data.srcRatio).abs > 1e-10
          ratio = @lastRatio + @outGen * (data.srcRatio - @lastRatio) / @outCount
        end

        floatIncrement = @indexInc * (ratio < 1.0 ? ratio : 1.0)
        increment = toFp(floatIncrement)

        startFilterIndex = toFp(inputIndex * floatIncrement)

        calcOutput(increment, startFilterIndex, floatIncrement / @indexInc, data.dataOut[@outGen..])
        @outGen += 2

        # Figure out the next index.
        inputIndex += (1.0 / ratio)
        rem = fmodOne(inputIndex)

        @bCurrent = ((@bCurrent.to_i64! + @channels * (inputIndex - rem).round.to_i64!) % @bLen).to_i32!
        inputIndex = rem
      end

      @lastPosition = inputIndex

      # Save current ratio rather then target ratio.
      @lastRatio = ratio

      data.inputFramesUsed = @inUsed.tdiv(@channels)
      data.outputFramesGen = @outGen.tdiv(@channels)
    end
  end

  # # A callback version of the `SincResamplerStereo` class.
  class SincResamplerStereoCb < SincResamplerStereo
    include CallbackResampler

    private def initialize(@channels : Int32)
      @coeffs = [] of Float32
      @buffer = [] of Float32
      @quality = Quality::Fast
      @callbackFunc = ->{ {Slice(Float32).empty, 0i64} }
      raise NotImplementedError.new("Use initialize(Int32, Float64, CallbackProc)")
    end

    # Creates a new `SincResamplerStereoCb` that will operate at the given
    # `Quality`.  The *ratio* parameter sets the initial ratio of the resampler
    # (`targetRate / sourceRate`), while *callbackFunc* is the method that will
    # be called whenever this resampler needs more input data.
    def initialize(@quality : Quality, ratio : Float64, @callbackFunc : CallbackProc)
      @channels = 2
      case @quality
      in .fast?
        @coeffs = LibRASincTables.sinc_fast_coeffs.to_unsafe
        @coeffHalfLen = LibRASincTables.sinc_fast_coeffs_size.to_i32 - 2
        @indexInc = LibRASincTables.sinc_fast_coeffs_index_inc.to_i32
      in .medium?
        @coeffs = LibRASincTables.sinc_medium_coeffs.to_unsafe
        @coeffHalfLen = LibRASincTables.sinc_medium_coeffs_size.to_i32 - 2
        @indexInc = LibRASincTables.sinc_medium_coeffs_index_inc.to_i32
      in .best?
        @coeffs = LibRASincTables.sinc_best_coeffs.to_unsafe
        @coeffHalfLen = LibRASincTables.sinc_best_coeffs_size.to_i32 - 2
        @indexInc = LibRASincTables.sinc_best_coeffs_index_inc.to_i32
      end

      @bLen = 3 * ((@coeffHalfLen + 2.0) / @indexInc * SRC_MAX_RATIO + 1).round.to_i32!
      @bLen = Math.max(@bLen, 4096)
      @bLen *= @channels
      @bLen += 1

      @buffer = Array(Float32).new(@bLen * @channels, 0.0f32)

      reset

      @data.srcRatio = ratio
      self.ratio = ratio
    end

    # This is not supported by this class since `LinearResamplerCb` uses the
    # `CallbackResampler` mixin.  Calling this always raises a
    # `NotImplementedError`.
    def process(source : Array(Float32)|Slice(Float32), dest : Array(Float32)|Slice(Float32),
                ratio : Float64) : Tuple(Int64, Int64, Bool)
      raise NotImplementedError.new("Not supported by callback resamplers")
    end

    # :ditto:
    def process(source : Array(Float64)|Slice(Float64), dest : Array(Float64)|Slice(Float64),
                ratio : Float64) : Tuple(Int64, Int64, Bool)
      raise NotImplementedError.new("Not supported by callback resamplers")
    end

    protected def processData : Nil
      variProcess(@data)
    end

    # :inherit:
    def reset : Nil
      super
      @data.reset
    end
  end
end