frequency_response.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/biquad.hpp>
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstddef>
#include <span>
 
namespace pulp::signal {

inline constexpr float min_response_db = -200.0f;

inline float magnitude_to_db(double magnitude) {
    if (!(magnitude > 0.0)) return min_response_db;
    double db = 20.0 * std::log10(magnitude);
    if (!std::isfinite(db) || db < min_response_db) return min_response_db;
    return static_cast<float>(db);
}

template <typename SampleType>
inline double section_magnitude(const BiquadCoefficientsT<SampleType>& c, double omega) {
    const std::complex<double> z_inv = std::polar(1.0, -omega);
    const std::complex<double> z_inv2 = z_inv * z_inv;
    const std::complex<double> num =
        static_cast<double>(c.b0) + static_cast<double>(c.b1) * z_inv + static_cast<double>(c.b2) * z_inv2;
    const std::complex<double> den =
        1.0 + static_cast<double>(c.a1) * z_inv + static_cast<double>(c.a2) * z_inv2;
    const double d = std::abs(den);
    if (!(d > 0.0)) return 0.0; // pole on the unit circle — treat as no output
    return std::abs(num) / d;
}
 
// Preserve the original concrete float entry point so aggregate calls such as
// section_magnitude({1, 0, 0, 0, 0}, omega) remain source-compatible.
inline double section_magnitude(const BiquadCoefficients& c, double omega) {
    return section_magnitude<float>(c, omega);
}

template <typename SampleType>
inline double section_group_delay(const BiquadCoefficientsT<SampleType>& c, double omega) {
    const std::complex<double> z1(std::cos(omega), -std::sin(omega));
    const std::complex<double> z2 = z1 * z1;
    const auto delay = [&](double p0, double p1, double p2) {
        return std::real((p1 * z1 + 2.0 * p2 * z2) / (p0 + p1 * z1 + p2 * z2));
    };
    return delay(c.b0, c.b1, c.b2) - delay(1.0, c.a1, c.a2);
}

template <typename SampleType, std::size_t Extent>
inline double cascade_magnitude(std::span<const BiquadCoefficientsT<SampleType>, Extent> sos,
                                double omega) {
    double mag = 1.0;
    for (const auto& c : sos) mag *= section_magnitude(c, omega);
    return mag;
}
 
// Preserve the original float-span entry point so containers such as
// std::vector<BiquadCoefficients> retain their implicit span conversion. The
// templated overload above adds precision without narrowing that source API.
inline double cascade_magnitude(std::span<const BiquadCoefficients> sos, double omega) {
    return cascade_magnitude<float, std::dynamic_extent>(sos, omega);
}

inline double angular_frequency(double freq_hz, double sample_rate) {
    if (!(sample_rate > 0.0)) return 0.0;
    const double omega = 2.0 * std::acos(-1.0) * freq_hz / sample_rate;
    return std::clamp(omega, 0.0, std::acos(-1.0));
}

template <typename SampleType>
inline float section_magnitude_db(const BiquadCoefficientsT<SampleType>& c, double freq_hz,
                                  double sample_rate) {
    return magnitude_to_db(section_magnitude(c, angular_frequency(freq_hz, sample_rate)));
}
 
inline float section_magnitude_db(const BiquadCoefficients& c, double freq_hz,
                                  double sample_rate) {
    return section_magnitude_db<float>(c, freq_hz, sample_rate);
}

template <typename SampleType, std::size_t Extent>
inline float cascade_magnitude_db(std::span<const BiquadCoefficientsT<SampleType>, Extent> sos,
                                  double freq_hz, double sample_rate) {
    return magnitude_to_db(cascade_magnitude(sos, angular_frequency(freq_hz, sample_rate)));
}
 
inline float cascade_magnitude_db(std::span<const BiquadCoefficients> sos, double freq_hz,
                                  double sample_rate) {
    return cascade_magnitude_db<float, std::dynamic_extent>(sos, freq_hz, sample_rate);
}

template <typename SampleType>
inline float magnitude_db(const BiquadT<SampleType>& filter, double freq_hz, double sample_rate) {
    const auto c = filter.coefficients();
    return section_magnitude_db(c, freq_hz, sample_rate);
}

inline double fft_bin_frequency(std::size_t index, std::size_t bin_count, double sample_rate) {
    if (bin_count <= 1) return 0.0;
    const double nyquist = sample_rate * 0.5;
    return nyquist * static_cast<double>(index) / static_cast<double>(bin_count - 1);
}

inline double log_frequency_at(std::size_t index, std::size_t count, double min_hz, double max_hz) {
    if (count <= 1) return min_hz;
    const double log_min = std::log10(min_hz);
    const double log_max = std::log10(max_hz);
    const double t = static_cast<double>(index) / static_cast<double>(count - 1);
    return std::pow(10.0, log_min + t * (log_max - log_min));
}

template <typename SampleType, std::size_t Extent>
inline void response_curve_db(std::span<const BiquadCoefficientsT<SampleType>, Extent> sos,
                              double min_hz, double max_hz, double sample_rate,
                              std::span<float> out) {
    const std::size_t n = out.size();
    for (std::size_t i = 0; i < n; ++i) {
        const double hz = log_frequency_at(i, n, min_hz, max_hz);
        out[i] = cascade_magnitude_db(sos, hz, sample_rate);
    }
}
 
inline void response_curve_db(std::span<const BiquadCoefficients> sos, double min_hz,
                              double max_hz, double sample_rate, std::span<float> out) {
    response_curve_db<float, std::dynamic_extent>(sos, min_hz, max_hz, sample_rate, out);
}
 
} // namespace pulp::signal