nth-order Recursive Infinite Impulse Response Filter - realtime-safe 64bit internal processing, change coeffs (and filter order) at runtime

filter-nth-order.cpp

/*
 * 
 *
 *   - `nth_order_system` Class: Manages the nth-order filter, including its
 * internal state and coefficients.
 *   - `compute_sample`: Computes a single output sample based on the input and
 * current state using Direct Form II Transposed structure.
 *   - `set_coefficients`: Sets new filter coefficients.
 *   - `reset`: Resets the internal state variables.
 *
 * - Real-time Configuration:
 *   - `handle_interactive_coefficient_changes`: Allows live updating of filter
 * coefficients via user input.
 *
 *
 * clang++ -std=c++20 filter-jack.cpp -o filter-jack -ljack -O2
 */

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <csignal>
#include <iostream>
#include <jack/jack.h>
#include <memory>
#include <numeric>
#include <sstream>
#include <string>

constexpr std::size_t max_order = 20;
using sample_type = double;class nth_order_system {
public:
  nth_order_system() { reset(); }

    //  single sample output 
    sample_type compute_sample(sample_type in) noexcept {
    double out;

    // to avoid repeated atomic loads
    auto a_local = a.load(std::memory_order_acquire);
    auto b_local = b.load(std::memory_order_acquire);

    double w_0 = in - std::inner_product(a_local->begin() + 1, a_local->begin() + a_size, z.begin(), 0.0);

    out = b_local->at(0) * w_0 + std::inner_product(b_local->begin() + 1, b_local->begin() + b_size, z.begin(), 0.0);

    std::rotate(z.rbegin(), z.rbegin() + 1, z.rend());
    z[0] = w_0;

    // Optional: Soft clipping 
    out = std::tanh(out);

    return static_cast<sample_type>(out);
  }

  void reset() noexcept { z.fill(0); }

  void set_coefficients(const std::array<double, max_order> &a_coeffs,
                        std::size_t a_size_,
                        const std::array<double, max_order + 1> &b_coeffs,
                        std::size_t b_size_) noexcept {

    auto new_a = std::make_shared<std::array<double, max_order>>(a_coeffs);
    auto new_b = std::make_shared<std::array<double, max_order + 1>>(b_coeffs);

    a.store(new_a, std::memory_order_release);
    b.store(new_b, std::memory_order_release);

    a_size = a_size_;
    b_size = b_size_;
    z_size = std::max(a_size, b_size) - 1;
    reset();
  }

private:
  std::atomic<std::shared_ptr<std::array<double, max_order>>> a{
      std::make_shared<std::array<double, max_order>>()};
  std::atomic<std::shared_ptr<std::array<double, max_order + 1>>> b{
      std::make_shared<std::array<double, max_order + 1>>()};
  std::array<double, max_order> z{}; // State variable
  std::size_t a_size{0};             // Size of 'a' coefficients
  std::size_t b_size{0};             // Size of 'b' coefficients
  std::size_t z_size{0};             // Required size of state variable 'z'
};


nth_order_system filter; // Global filter instance

// JACK Audio process callback function
int process(jack_nframes_t nframes, void *arg) {
  auto *ports = static_cast<std::pair<jack_port_t *, jack_port_t *> *>(arg);
  auto *in = static_cast<jack_default_audio_sample_t *>(
      jack_port_get_buffer(ports->first, nframes));
  auto *out = static_cast<jack_default_audio_sample_t *>(
      jack_port_get_buffer(ports->second, nframes));

  for (unsigned int i = 0; i < nframes; ++i) {
    double computed_sample = filter.compute_sample(static_cast<double>(in[i]));
    out[i] = static_cast<jack_default_audio_sample_t>(computed_sample);
  }

  return 0;
}

jack_client_t *initialize_jack(const char *client_name) {
  jack_options_t options = JackNullOption;
  jack_status_t status;
  jack_client_t *client = jack_client_open(client_name, options, &status);

  if (!client) {
    std::cerr << "jack_client_open() failed, status = 0x" << std::hex << status
              << '\n';
    if (status & JackServerFailed) {
      std::cerr << "Unable to connect to JACK server\n";
    }
    return nullptr;
  }

  if (status & JackNameNotUnique) {
    client_name = jack_get_client_name(client);
    std::cerr << "Unique name assigned: " << client_name << '\n';
  }

  return client;
}

void configure_filter() {
  std::array<double, max_order> a_coeffs = {0.2, 0.4, 0.3, 0.1};
  std::array<double, max_order + 1> b_coeffs = {1.0, 0.7, 0.5, 0.3, 0.1};
  filter.set_coefficients(a_coeffs, 4, b_coeffs, 5);
}

// Register input and output ports with JACK
std::pair<jack_port_t *, jack_port_t *> register_ports(jack_client_t *client) {
  jack_port_t *input_port = jack_port_register(
      client, "input", JACK_DEFAULT_AUDIO_TYPE, JackPortIsInput, 0);
  if (!input_port) {
    std::cerr << "No more JACK ports available\n";
    return {nullptr, nullptr};
  }

  jack_port_t *output_port = jack_port_register(
      client, "output", JACK_DEFAULT_AUDIO_TYPE, JackPortIsOutput, 0);
  if (!output_port) {
    std::cerr << "No more JACK ports available\n";
    return {nullptr, nullptr};
  }

  return {input_port, output_port};
}

bool set_process_callback(jack_client_t *client,
                          std::pair<jack_port_t *, jack_port_t *> &ports) {
  if (jack_set_process_callback(client, process, &ports) != 0) {
    std::cerr << "jack_set_process_callback() failed\n";
    return false;
  }
  return true;
}

bool activate_client(jack_client_t *client) {
  if (jack_activate(client) != 0) {
    std::cerr << "Cannot activate client\n";
    return false;
  }
  return true;
}

std::atomic<bool> running(true); // Atomic flag to handle termination signal

void signal_handler(int signal) {
  (void)signal;
  running.store(false);
}

void handle_interactive_coefficient_changes() {
  std::cout
      << "Instructions: Enter new coefficients as space-separated values.\n";
  std::cout
      << "First, enter the 'a' coefficients, then the 'b' coefficients.\n";
  std::cout << "Type 'quit' to exit the program.\n";

  std::string input;

  while (running.load()) {
    std::array<double, max_order> new_a_coeffs{};
    std::array<double, max_order + 1> new_b_coeffs{};
    std::size_t a_count = 0, b_count = 0;

    std::cout << "Enter 'a' coefficients (a0 a1 ... an) or 'quit' to exit: ";
    std::getline(std::cin, input);
    if (input == "quit") {
      running.store(false);
      break;
    }

    std::istringstream iss_a(input);
    double coeff;
    bool valid_input = true;

    while (iss_a >> coeff) {
      if (a_count < max_order) {
        new_a_coeffs[a_count++] = coeff;
      } else {
        std::cerr << "Too many 'a' coefficients entered. Maximum allowed is "
                  << max_order << ".\n";
        valid_input = false;
        break;
      }
    }

    if (!valid_input || a_count == 0) {
      continue;
    }

    std::cout << "Enter 'b' coefficients (b0 b1 ... bn) or 'quit' to exit: ";
    std::getline(std::cin, input);
    if (input == "quit") {
      running.store(false);
      break;
    }

    std::istringstream iss_b(input);

    valid_input = true;
    while (iss_b >> coeff) {
      if (b_count < max_order + 1) {
        new_b_coeffs[b_count++] = coeff;
      } else {
        std::cerr << "Too many 'b' coefficients entered. Maximum allowed is "
                  << (max_order + 1) << ".\n";
        valid_input = false;
        break;
      }
    }

    if (!valid_input || b_count == 0) {
      continue;
    }

    try {
      filter.set_coefficients(new_a_coeffs, a_count, new_b_coeffs, b_count);
      std::cout << "Coefficients updated.\n";
    } catch (const std::exception &e) {
      std::cerr << "Failed to update coefficients: " << e.what() << '\n';
    }
  }
}

int main(int argc, char *argv[]) {
  const char *client_name = "nth_order_filter";

  if (argc > 1) {
    client_name = argv[1];
  }

  jack_client_t *client = initialize_jack(client_name);
  if (!client)
    return 1;

  configure_filter();

  auto ports = register_ports(client);
  if (!ports.first || !ports.second)
    return 1;

  if (!set_process_callback(client, ports))
    return 1;

  if (!activate_client(client))
    return 1;

  std::signal(SIGINT, signal_handler);

  handle_interactive_coefficient_changes();

  jack_client_close(client);
  return 0;
}