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

###
### https://www.musicdsp.org/en/latest/Filters/197-rbj-audio-eq-cookbook.html
### plus some referencing of the SoX source code (biquads.c)
###

module RemiAudio::DSP
  # A digital biquad lowpass filter.
  #
  # This is identical to using a normal `BiQuadFilter` instance, except it is
  # optimized strictly for lowpass usage and includes the `LPFilter` module.
  class BiquadLP
    include LPFilter

    # :nodoc:
    RESONANCE_PEAK_OFFSET = 1.0f64 - (1.0f64 / Math.sqrt(2.0f64))

    getter b0 : Float64 = 0.0f64
    getter b1 : Float64 = 0.0f64
    getter b2 : Float64 = 0.0f64
    getter a1 : Float64 = 0.0f64
    getter a2 : Float64 = 0.0f64

    # The current frequency of the filter.
    @[AlwaysInline]
    getter cutoff : Float64 = CUTOFF_MAX

    # The current Q value of the filter.
    @[AlwaysInline]
    getter resonance : Float64 = Float64::MIN_POSITIVE

    @x1 : Float64 = 0.0f64
    @x2 : Float64 = 0.0f64
    @y1 : Float64 = 0.0f64
    @y2 : Float64 = 0.0f64

    @sampleRate : Float64
    @invSampleRate : Float64

    # Creates a new `BiQuadFilter`.
    def initialize(newSampleRate)
      @sampleRate = newSampleRate.to_f64!
      @invSampleRate = 1.0 / @sampleRate
      set(CUTOFF_MAX, RESONANCE_MIN)
    end

    # Clears the internal buffer.  This does not change the settings of the EQ.
    @[AlwaysInline]
    def reset : Nil
      @x1 = 0.0f64
      @x2 = 0.0f64
      @y1 = 0.0f64
      @y2 = 0.0f64
    end

    # :nodoc:
    RES_ADJUST = RemiAudio.decibelsToLinear(48.0)

    @[AlwaysInline]
    def updateCoefficients : Nil
      omega : Float64 = RemiMath::TWO_PI * @cutoff * @invSampleRate
      @active = @cutoff < 0.499 * @sampleRate # Was <= ENABLE_AT
      res = RemiAudio.decibelsToLinear(@resonance * RES_ADJUST)
      q = res - RESONANCE_PEAK_OFFSET / (1.0 + 6.0 * (res - 1.0))
      alpha : Float64 = Math.sin(omega) / (2.0 * q)
      cosO  : Float64 = Math.cos(omega)
      @b0 = (1.0f64 - cosO) * 0.5f64
      @b1 = 1.0f64 - cosO
      @b2 = (1.0f64 - cosO) * 0.5f64
      @a1 = -2.0f64 * cosO
      @a2 = 1.0f64 - alpha

      a0 : Float64 = 1.0f64 + alpha
      @b0 /= a0
      @b1 /= a0
      @b2 /= a0
      @a1 /= a0
      @a2 /= a0
    end

    def cutoff=(val : Float64) : Nil
      return if val == @cutoff
      @cutoff = val.clamp(CUTOFF_MIN, CUTOFF_MAX)
      updateCoefficients
    end

    def resonance=(val : Float64) : Nil
      return if @resonance == val
      @resonance = val.clamp(Float64::MIN_POSITIVE, RESONANCE_MAX)
      updateCoefficients
    end

    def set(newCutoff : Float64, newResonance : Float64)
      return if @cutoff == newCutoff && @resonance == newResonance
      @cutoff = newCutoff.clamp(CUTOFF_MIN, CUTOFF_MAX)
      @resonance = newResonance.clamp(Float64::MIN_POSITIVE, RESONANCE_MAX - 0.99)
      updateCoefficients
    end

    @[AlwaysInline]
    def process(sample : Float64) : Float64
      return sample unless @active
      output : Float64 = ((@b0 * sample) + (@b1 * @x1) + (@b2 * @x2)) - (@a1 * @y1) - (@a2 * @y2)
      @x2 = @x1
      @x1 = sample
      @y2 = @y1
      @y1 = output
    end

    # Generates a string that can be passed to a program to plot this filter on
    # a graph.  The `PlotType` dictates what kind of script is generated.
    def plot(pt : PlotType, *, coeffsOnly : Bool = false) : String
      # Adapted from what SoX prints out.
      String.build do |str|
        case pt
        in .gnu_plot?
          if coeffsOnly
            str << sprintf("b0 = %.15e; ", @b0)
            str << sprintf("b1 = %.15e; ", @b1)
            str << sprintf("b2 = %.15e; ", @b2)
            str << sprintf("a1 = %.15e; ", @a1)
            str << sprintf("a2 = %.15e", @a2)
          else
            str << "# gnuplot file\n"
            str << "set title 'BiQuad Filter, sample rate: #{@sampleRate}'\n"
            str << "set xlabel 'Frequency (Hz)'\n"
            str << "set ylabel 'Amplitude Response (dB)'\n"
            str << sprintf("Fs = %g\n", @sampleRate)
            str << sprintf("b0 = %.15e; ", @b0)
            str << sprintf("b1 = %.15e; ", @b1)
            str << sprintf("b2 = %.15e; ", @b2)
            str << sprintf("a1 = %.15e; ", @a1)
            str << sprintf("a2 = %.15e\n", @a2)
            str << "o = 2 * pi / Fs\n"
            str << "H(f) = sqrt((b0 * b0 + b1 * b1 + b2 * b2 + 2.0 * (b0 * b1 + b1 * b2) * " <<
              "cos(f * o) + 2.0 * (b0 * b2) * cos(2.0 * f * o)) / " <<
              "(1.0 + a1 * a1 + a2 * a2 + 2.0 * (a1 + a1 * a2) * " <<
              "cos(f * o) + 2.0 * a2 * cos(2.0 * f * o)))\n"
            str << "set logscale x\n"
            str << "set samples 250\n"
            str << "set grid xtics ytics\n"
            str << "set key off\n"
            str << "plot [f=10 : Fs/2] [-35:25] 20 * log10(H(f))"
          end
        in .octave?
          str << "%% GNU Octave file (may also work with MATLAB(R) )\n"
          str << "Fs = #{@sampleRate};\n"
          str << "minF = 10;\n"
          str << "maxF = Fs/2;\n"
          str << "sweepF = logspace(log10(minF), log10(maxF), 200);\n"
          str << sprintf("[h, w] = freqz([%.15e %.15e %.15e], [1 %.15e %.15e], sweepF, Fs);\n",
                         @b0, @b1, @b2, @a1, @a2)
          str << "semilogx(w, 20 * log10(h))\n"
          str << "title('BiQuad Filter, sample rate: #{@sampleRate}')\n"
          str << "xlabel('Frequency (Hz)')\n"
          str << "ylabel('Amplitude Response (dB)')\n"
          str << "axis([minF maxF -35 25])\n"
          str << "grid on\n"
          str << "disp('Hit return to continue')\n"
          str << "pause\n"
        end
      end
    end
  end
end
