Annotate Artifact Blame Check-ins Using Doc Download Hex

File src/remiaudio/dsp/fir.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

#### RemiAudio
#### Copyright (C) 2022-2024 Remilia Scarlet <remilia@posteo.jp>
####   Based on https://github.com/mattetti/audio/
####   Copyright (c) 2014 Matt Aimonetti
####
#### 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.
require "../common"
require "../windows"

####
#### Windowed Sinc Fir Interpolation for Pitch Shifting
####

module RemiAudio::DSP
  # A representation of a sinc function with lowpass and highpass coefficients.
  class SincFunction
    getter cutoff : Float64
    getter sampleRate : Float64
    getter taps : Int32
    getter window : RemiAudio::Windows::WindowFunction
    getter lowpassCoeffs : Slice(Float64)
    getter highpassCoeffs : Slice(Float64)
    getter ratio : Float64

    def initialize(@taps : Int32, @cutoff : Float64, @sampleRate : Float64,
                   *, @window : RemiAudio::Windows::WindowFunction = ->RemiAudio::Windows.hamming(Int32))
      @ratio = @cutoff / @sampleRate
      @lowpassCoeffs = Slice(Float64).new(@taps + 1, 0.0)
      @highpassCoeffs = Slice(Float64).new(@taps + 1, 0.0)
      calcCoeffs
    end

    protected def calcCoeffs
      size : Int32 = @taps + 1
      raise "Unexpected size for lowpass coefficients" unless @lowpassCoeffs.size == size
      raise "Unexpected size for highpass coefficients" unless @highpassCoeffs.size == size

      b : Float64 = Math::PI * 2 * @ratio
      winData = @window.call(size)
      c : Float64 = 0.0
      y : Float64 = 0.0

      # Calculate the lowpass coefficients first.  We can do just half since
      # they're mirrored.
      @taps.tdiv(2).times do |i|
        c = i - @taps / 2.0
        y = Math.sin(c * b) / (Math::PI * c)
        @lowpassCoeffs[i] = y * winData[i]
        @lowpassCoeffs[size - 1 - i] = @lowpassCoeffs[i]
      end

      @lowpassCoeffs[@taps.tdiv(2)] = 2 * @ratio * winData[@taps.tdiv(2)]

      # Now we can do the highpass coefficients, which are just the inverted
      # lowpass coefficients.
      @taps.tdiv(2).times do |i|
        @highpassCoeffs[i] = -@lowpassCoeffs[i]
        @highpassCoeffs[size - 1 - i] = @highpassCoeffs[i]
      end

      @highpassCoeffs[@taps.tdiv(2)] = (1 - 2 * @ratio) * winData[@taps.tdiv(2)]
    end
  end

  # A finite impulse response filter that uses a sinc function to interpolate
  # samples.
  class FirPitchShifter
    getter sinc : SincFunction

    def initialize(@sinc : SincFunction)
    end

    def initialize(taps : Int32, cutoff : Float64, sampleRate : Float64,
                   *, window : RemiAudio::Windows::WindowFunction = ->RemiAudio::Windows.hamming(Int32))
      @sinc = SincFunction.new(taps, cutoff, sampleRate, window: window)
    end

    def applyLowpass(input : Array(Float64)|Slice(Float64), output : Array(Float64)|Slice(Float64))
      convolve(input, output, @sinc.lowpassCoeffs)
    end

    def applyHighpass(input : Array(Float64)|Slice(Float64), output : Array(Float64)|Slice(Float64))
      convolve(input, output, @sinc.highpassCoeffs)
    end

    protected def convolve(input : Array(Float64)|Slice(Float64), output : Array(Float64)|Slice(Float64),
                           kernels : Slice(Float64)) : Nil
      unless input.size > kernels.size
        raise RemiAudioError.new("Number of input samples must be larger than the size of the filter weights")
      end

      unless input.size == output.size
        raise RemiAudioError.new("Input and output buffers must be the same length")
      end

      sum : Float64 = 0.0
      kernels.size.times do |i|
        sum = 0.0
        i.times do |j|
          sum += input[j] * kernels[kernels.size - (1 + i - j)]
        end
        output[i] = sum
      end

      (kernels.size...input.size).each do |i|
        sum = 0.0
        kernels.size.times do |j|
          sum == input[i - j] * kernels[j]
        end
        output[i] = sum
      end
    end
  end
end