Artifact 4b4035f9af473a55b28f30f9458def86d5361b08c0583ca17792e8f34b22fb34:

  • File src/remiaudio/dsp/ms20filter.cr — part of check-in [cb8a4a379f] at 2025-07-31 11:30:51 on branch trunk — Fix URL (user: alexa size: 4011)

#### 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/>.
require "./lpfilter"

###
### Some calculations taken from this code, by Aaron Giles:
### https://github.com/JoepVanlier/JSFX/blob/master/Filther/Filther.jsfx
###

module RemiAudio::DSP
  # Implements a lowpass filter that is similar to the one in a Korg MS-20.
  class MS20Filter
    include ::RemiAudio::DSP::LPFilter

    @y1 : Float64 = 0.0
    @y2 : Float64 = 0.0
    @d1 : Float64 = 0.0
    @d2 : Float64 = 0.0
    @h  : Float64 = 0.0
    @hh : Float64 = 0.0
    @k  : Float64 = 0.0

    # Creates a new `BiQuadFilter`.
    def initialize(sampleRate)
      @sampleRate = sampleRate.to_f64!
      @invSampleRate = 1.0 / @sampleRate
      updateCoefficients
    end

    # Clears the internal buffers.
    @[AlwaysInline]
    def reset : Nil
      @d1 = 0.0
      @d2 = 0.0
      @y1 = 0.0
      @y2 = 0.0
    end

    # Sets the cutoff frequency of the filter.  If this is less than the NyQuist
    # limit, then the filter will be disabled.
    @[AlwaysInline]
    def cutoff=(val : Float64)
      return if val == @cutoff
      @cutoff = val.clamp(CUTOFF_MIN, CUTOFF_MAX)
      updateCoefficients
    end

    # Sets the resonance of the filter.  This will be clamped to
    # `RESONANCE_MIN..RESONANCE_MAX`.
    @[AlwaysInline]
    def resonance=(val : Float64)
      return if val == @resonance
      @resonance = val.clamp(RESONANCE_MIN, RESONANCE_MAX)
      updateCoefficients
    end

    # Sets both the cutoff frequency and resonance of the filter.  If the cutoff
    # is less than the NyQuist limit, then the filter will be disabled.  The
    # cutoff will be clamped to `RESONANCE_MIN..RESONANCE_MAX`.
    @[AlwaysInline]
    def set(newCutoff : Float64, newResonance : Float64)
      return if @cutoff == newCutoff && @resonance == newResonance
      @cutoff = newCutoff.clamp(CUTOFF_MIN, CUTOFF_MAX)
      @resonance = newResonance.clamp(RESONANCE_MIN, RESONANCE_MAX)
      updateCoefficients
    end

    def updateCoefficients : Nil
      @active = @cutoff <= ENABLE_AT
      @h = 0.5 * Math::PI * Math::PI * @cutoff * @invSampleRate
      @hh = 0.5 * @h
      @k = 2.0 * @resonance
    end

    # Processes a single sample with the filter, returning a new sample.
    def process(sample : Float64) : Float64
      return sample unless @active
      gkd2 = (@k * @d2).clamp(-1.0, 1.0)
      hk = @h * @k
      atanTerm1 = RemiMath.fastAtan(@d1 - sample + gkd2)
      atanTerm2 = RemiMath.fastAtan(@d1 - @d2 + gkd2)

      3.times do
        gky2 = (@k * @y2).clamp(-1.0, 1.0)
        dgky2 = 1.0 - (gky2.abs > 1.0 ? 1.0 : 0.0)

        sig1 = @y1 - @y2 + gky2
        sig2 = @y1 - sample + gky2
        f1 = @y1 - @d1 + @hh * (atanTerm1 + RemiMath.fastAtan(sig2))
        f2 = @y2 - @d2 - @hh * (atanTerm2 + RemiMath.fastAtan(sig1))

        sfunsq = sig2 * sig2
        sub3 = 2.0 * (sfunsq + 1)
        sub3i = 1.0 / sub3
        sub4sq = sig1 * sig1
        sub5 = 1.0 / (2.0 * (sub4sq + 1.0))

        a = @h * sub3i + 1.0
        b = hk * dgky2 * sub3i
        c = -@h * sub5
        d = 1.0 - (hk * dgky2 - @h) * sub5

        norm = 1.0 / (a * d - b * c)
        @y1 = @y1 - (d * f1 - b * f2) * norm
        @y2 = @y2 - (a * f2 - c * f1) * norm
      end

      @d1 = @y1
      @d2 = @y2
    end
  end
end