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