LibMedia: Implement biquad filter frequency response

Implement a generic function for a biquad filter frequency
response for use in LibWeb/WebAudio. Also add the computation
for the low pass filter coefficients for testing.

AK/Complex.h

@@ -8,6 +8,7 @@
 
 #include <AK/Concepts.h>
 #include <AK/Math.h>
+#include <AK/String.h>
 
 namespace AK {
 

Tests/LibMedia/CMakeLists.txt

@@ -2,6 +2,7 @@ set(TEST_SOURCES
     TestH264Decode.cpp
     TestParseMatroska.cpp
     TestPlaybackStream.cpp
+    TestSignalProcessing.cpp
     TestVorbisDecode.cpp
     TestVP9Decode.cpp
     TestWav.cpp

Tests/LibMedia/TestSignalProcessing.cpp

@@ -0,0 +1,43 @@
+/*
+ * Copyright (c) 2024, Noah Bright <[email protected]>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#include <AK/Complex.h>
+#include <LibTest/TestCase.h>
+
+#include <LibMedia/Audio/SignalProcessing.h>
+
+TEST_CASE(test_biquad_filter_frequency_response)
+{
+    size_t n_frequencies = 4096;
+    Vector<f64> frequencies;
+    frequencies.resize(n_frequencies);
+
+    // roughly the range of human hearing in Hz
+    auto delta_f = (20000 - 20) / n_frequencies;
+    for (size_t i = 0; i < n_frequencies; i++)
+        frequencies[i] = 20.0 + (f64)i * (f64)delta_f;
+
+    // default coefficients from https://webaudio.github.io/web-audio-api/#BiquadFilterOptions
+    f64 Q = 1.0;
+    f64 frequency = 350;
+    f64 sample_rate = 44100;
+
+    auto omega_0 = 2 * AK::Pi<f64> * frequency / sample_rate;
+    auto alpha_Q_dB = sin(omega_0) / (2 * pow(10, Q / 20));
+
+    auto lowpass_coeffs = Audio::biquad_filter_lowpass_coefficients(omega_0, alpha_Q_dB);
+
+    auto frequency_response = Audio::biquad_filter_frequency_response(frequencies, lowpass_coeffs);
+    for (auto response : frequency_response) {
+        auto phase_response = response.phase();
+        VERIFY(phase_response != AK::NaN<f64>);
+        VERIFY(phase_response != AK::Infinity<f64>);
+
+        auto magnitude_response = response.magnitude();
+        VERIFY(magnitude_response != AK::NaN<f64>);
+        VERIFY(magnitude_response != AK::Infinity<f64>);
+    }
+}

Userland/Libraries/LibMedia/Audio/SignalProcessing.h

@@ -0,0 +1,70 @@
+/*
+ * Copyright (c) 2024, Noah Bright <[email protected]>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#pragma once
+
+#include <AK/Array.h>
+#include <AK/Complex.h>
+#include <AK/Vector.h>
+
+namespace Audio {
+
+template<AK::Concepts::Arithmetic T>
+Vector<Complex<T>> biquad_filter_frequency_response(Vector<T> const& frequencies, Array<T, 6> const& coefficients)
+{
+    // frequencies is expected to be frequencies in Hz, not angular frequencies in rad/sec
+    // coefficients should be a fixed sized array of 6 containing [b0, b1, b2, a0, a1, a2]
+    //
+    // the transfer function of a biquadratic filter is
+    //      H(z) = ( b0/a0 + b1/a0 * z^-1 + b2/a0 * z^-2 ) / ( 1 + a1/a0 * z^-1 + a2/a0 * z^-2)
+    //
+    // as written, the numerator and denominator above both have 3 floating point multiplications
+    // rewriting the numerator, that can be reduced to 2:
+    // b0/a0 + (b1/a0 + b2/a0 * z) * z^-1
+    //
+    // the frequency response of a filter at frequency omega is its transfer function evaluated
+    // at z = e^{i omega}, with angular frequency omega = 2*pi*f
+    //
+    // the magnitude(phase) response of a filter is the magnitude(phase) of its frequency response
+
+    T B0 = coefficients[0] / coefficients[3];
+    T B1 = coefficients[1] / coefficients[3];
+    T B2 = coefficients[2] / coefficients[3];
+    T A1 = coefficients[4] / coefficients[3];
+    T A2 = coefficients[5] / coefficients[3];
+
+    auto transfer_function = [&](Complex<T> z) {
+        auto z_inv = 1 / z;
+        auto numerator = B0 + (B1 + B2 * z_inv) * z_inv;
+        auto denominator = 1.0 + (A1 + A2 * z_inv) * z_inv;
+        return numerator / denominator;
+    };
+
+    Vector<Complex<T>> frequency_response;
+    frequency_response.ensure_capacity(frequencies.size());
+
+    for (size_t k = 0; k < frequency_response.size(); k++) {
+        auto i = Complex<T>(0, 1);
+        auto arg = cexp(2 * AK::Pi<T> * i * frequencies[k]);
+        frequency_response[k] = transfer_function(arg);
+    }
+
+    return frequency_response;
+}
+
+template<AK::Concepts::Arithmetic T>
+Array<T, 6> biquad_filter_lowpass_coefficients(T omega_0, T alpha_Q_dB)
+{
+    auto b0 = 0.5 * (1 - cos(omega_0));
+    auto b1 = 1 - cos(omega_0);
+    auto b2 = b0;
+    auto a0 = 1 + alpha_Q_dB;
+    auto a1 = -2 * cos(omega_0);
+    auto a2 = 1 - alpha_Q_dB;
+    return { b0, b1, b2, a0, a1, a2 };
+}
+
+}