[PERFORMANCE] Optimize long linear convolution scheduling

NAM/linear.h

@@ -7,6 +7,12 @@ namespace nam
 
 struct LinearFFTState;
 
+struct LinearFFTPlan
+{
+  int direct_taps;
+  int max_partition_size;
+};
+
 /// \brief Selects the convolution engine used by Linear models.
 enum class LinearImplementation
 {
@@ -62,7 +68,10 @@ class Linear : public Buffer
   void _configure_fft_state();
   void _process_direct(NAM_SAMPLE** input, NAM_SAMPLE** output, const int num_frames);
   void _process_fft(NAM_SAMPLE** input, NAM_SAMPLE** output, const int num_frames);
-  void _run_fft_block(const int channel);
+  void _advance_fft_job(const int tier, const int channel);
+  void _advance_fft_jobs(const int channel);
+  void _start_fft_block(const int tier, const int channel, const long long block_start);
+  void _finish_fft_block(const int tier, const int channel);
 };
 
 namespace linear
@@ -86,6 +95,12 @@ LinearImplementation parse_implementation(const std::string& implementation);
 /// \brief String name for a Linear implementation.
 std::string implementation_to_string(const LinearImplementation implementation);
 
+/// \brief Select the tuned convolution plan for an impulse-response length.
+LinearFFTPlan select_fft_plan(int receptive_field);
+
+/// \brief Select the default implementation for an impulse-response length.
+LinearImplementation select_implementation(int receptive_field);
+
 /// \brief Parse Linear configuration from JSON
 /// \param config JSON configuration object
 /// \return LinearConfig

tools/BENCHMARK_LINEAR.md

@@ -0,0 +1,52 @@
+# Linear convolution benchmark
+
+`bench_linear` measures callback-time distribution for synthetic linear models. It intentionally reports callback
+times rather than only throughput: a convolution implementation can have a low average cost and still cause audio
+dropouts when periodic work exceeds the callback deadline.
+
+Build and run a Release benchmark with:
+
+```sh
+cmake -S . -B build-release -DCMAKE_BUILD_TYPE=Release
+cmake --build build-release --target bench_linear -j
+build-release/tools/bench_linear
+```
+
+The optional arguments are `taps`, `callback size`, `seconds`, `fft|direct`, and `cold|verify`. The `cold` mode
+touches a 64 MiB buffer before each timed callback to expose cache-sensitive plans. `verify` renders an impulse
+through the entire requested filter—including a one-minute filter—and compares every output sample with its tap.
+
+## Dispatch tuning
+
+The dispatch table was tuned on an Apple M1 at 48 kHz using the same `-O3` optimization level as the Release plugin.
+The main constraint was the 666.67 microsecond deadline of a 32-sample callback. Candidate plans varied direct-head
+sizes from 64 through 256 samples and maximum FFT partitions from 1,024 through 32,768 samples.
+
+The selected FFT plans produced these representative 10-second results (the default dispatcher uses direct
+convolution for the 1,024-tap case):
+
+| Taps | Direct head | Maximum partition | Mean (us) | p99 (us) | Max (us) | Overruns |
+|---:|---:|---:|---:|---:|---:|---:|
+| 1,024 | 128 | 256 | 2.37 | 10.62 | 55.67 | 0 |
+| 2,048 | 128 | 512 | 2.90 | 14.54 | 63.38 | 0 |
+| 4,096 | 128 | 1,024 | 3.65 | 23.79 | 58.33 | 0 |
+| 8,192 | 128 | 2,048 | 4.63 | 49.12 | 81.79 | 0 |
+| 48,000 | 64 | 4,096 | 6.45 | 48.83 | 123.42 | 0 |
+| 240,000 | 64 | 8,192 | 8.21 | 91.12 | 218.33 | 0 |
+| 1,200,000 | 64 | 8,192 | 14.41 | 99.79 | 297.58 | 0 |
+| 2,880,000 | 64 | 8,192 | 24.22 | 112.00 | 324.62 | 0 |
+
+The old 1,024-sample uniform implementation measured about 2.7--3.1 ms at p99 for the 1,208,121-tap atmospheric
+model and overran 32- and 64-sample callback deadlines. The non-uniform plan reduces that periodic burst by more than
+an order of magnitude.
+
+The table favors bounded callback time over the absolute minimum mean:
+
+- A 64-sample direct head had lower p99 times for long filters; a 128-sample head reduced overhead for 2K--8K filters.
+- A 16,384-sample maximum partition slightly reduced mean time for multi-second filters but roughly doubled p99 and
+  produced deadline outliers. A 4,096 maximum kept transforms smaller but nearly doubled one-minute steady-state cost.
+- Direct convolution through 1,024 taps had a flatter callback profile than FFT convolution at essentially the same
+  mean cost, so `Auto` selects it for that range.
+
+These numbers are hardware-specific. When changing the FFT backend, data layout, or dispatch table, rerun warm and
+cold-cache tests on every supported architecture and optimize for the worst callback distribution, not only mean CPU.

tools/CMakeLists.txt

@@ -18,6 +18,18 @@ include_directories(tools ${AUDIO_DSP_TOOLS_DIR}/dsp)
 
 add_executable(loadmodel loadmodel.cpp ${NAM_SOURCES})
 add_executable(benchmodel benchmodel.cpp ${NAM_SOURCES})
+add_executable(bench_linear bench_linear.cpp ${NAM_SOURCES})
+target_compile_features(bench_linear PUBLIC cxx_std_20)
+set_target_properties(bench_linear PROPERTIES
+	CXX_VISIBILITY_PRESET hidden
+	INTERPROCEDURAL_OPTIMIZATION TRUE
+	PREFIX ""
+)
+if (MSVC)
+	target_compile_options(bench_linear PRIVATE "$<$<CONFIG:RELEASE>:/O2>")
+else()
+	target_compile_options(bench_linear PRIVATE "$<$<CONFIG:RELEASE>:-O3>")
+endif()
 add_executable(render render.cpp ${NAM_SOURCES} ${AUDIO_DSP_TOOLS_WAV_SOURCES})
 target_compile_features(render PUBLIC cxx_std_20)
 # AudioDSPTools wav.cpp has sign-compare issues; don't fail build

tools/bench_linear.cpp

@@ -0,0 +1,128 @@
+#include "NAM/linear.h"
+
+#include <algorithm>
+#include <chrono>
+#include <cmath>
+#include <cstdlib>
+#include <iomanip>
+#include <iostream>
+#include <numeric>
+#include <string>
+#include <vector>
+
+namespace
+{
+using Clock = std::chrono::steady_clock;
+
+std::vector<float> make_weights(const int taps)
+{
+  std::vector<float> weights(taps);
+  for (int i = 0; i < taps; ++i)
+    weights[i] = 0.001f * std::exp(-6.0f * (float)i / (float)taps) * std::sin(0.031f * (float)(i + 1));
+  return weights;
+}
+
+double percentile(const std::vector<double>& sorted, const double fraction)
+{
+  const size_t index = std::min(sorted.size() - 1, (size_t)std::floor(fraction * (double)(sorted.size() - 1)));
+  return sorted[index];
+}
+
+bool verify_impulse_response(const int taps, const int callback_size, const int sample_rate)
+{
+  auto weights = make_weights(taps);
+  nam::Linear model(1, 1, taps, false, weights, sample_rate, nam::LinearImplementation::FFT);
+  model.Reset(sample_rate, callback_size);
+
+  std::vector<NAM_SAMPLE> input(callback_size, 0.0);
+  std::vector<NAM_SAMPLE> output(callback_size, 0.0);
+  NAM_SAMPLE* inputs[] = {input.data()};
+  NAM_SAMPLE* outputs[] = {output.data()};
+  double max_abs_error = 0.0;
+  int sample = 0;
+  const int samples_to_process = taps + 2 * callback_size;
+  while (sample < samples_to_process)
+  {
+    std::fill(input.begin(), input.end(), 0.0);
+    if (sample == 0)
+      input[0] = 1.0;
+    const int frames = std::min(callback_size, samples_to_process - sample);
+    model.process(inputs, outputs, frames);
+    for (int i = 0; i < frames; ++i)
+    {
+      const double expected = sample + i < taps ? weights[sample + i] : 0.0;
+      max_abs_error = std::max(max_abs_error, std::abs((double)output[i] - expected));
+    }
+    sample += frames;
+  }
+  std::cout << "verified_taps=" << taps << ",max_abs_error=" << std::scientific << max_abs_error << '\n';
+  return max_abs_error < 5.0e-5;
+}
+
+void run(const int taps, const int callback_size, const int sample_rate, const double seconds,
+         const nam::LinearImplementation implementation, const bool cold_cache)
+{
+  constexpr int max_callback_size = 512;
+  auto weights = make_weights(taps);
+  nam::Linear model(1, 1, taps, false, weights, sample_rate, implementation);
+  model.Reset(sample_rate, max_callback_size);
+
+  std::vector<NAM_SAMPLE> input(max_callback_size);
+  std::vector<NAM_SAMPLE> output(max_callback_size);
+  for (int i = 0; i < max_callback_size; ++i)
+    input[i] = (NAM_SAMPLE)(0.1 * std::sin(0.017 * i) + 0.03 * std::cos(0.043 * i));
+  NAM_SAMPLE* inputs[] = {input.data()};
+  NAM_SAMPLE* outputs[] = {output.data()};
+
+  for (int i = 0; i < 4096 / callback_size; ++i)
+    model.process(inputs, outputs, callback_size);
+
+  const int callbacks = std::max(1, (int)std::ceil(seconds * sample_rate / callback_size));
+  std::vector<double> durations;
+  durations.reserve(callbacks);
+  std::vector<unsigned char> cache_polluter(cold_cache ? 64 * 1024 * 1024 : 0, 1);
+  volatile NAM_SAMPLE sink = 0.0;
+  volatile unsigned int cache_sink = 0;
+  for (int i = 0; i < callbacks; ++i)
+  {
+    for (size_t byte = 0; byte < cache_polluter.size(); byte += 64)
+      cache_sink = cache_sink + cache_polluter[byte];
+    const auto start = Clock::now();
+    model.process(inputs, outputs, callback_size);
+    const auto end = Clock::now();
+    sink = sink + output[i % callback_size];
+    durations.push_back(std::chrono::duration<double, std::micro>(end - start).count());
+  }
+
+  std::sort(durations.begin(), durations.end());
+  const double deadline = 1.0e6 * callback_size / sample_rate;
+  const auto overruns =
+    std::count_if(durations.begin(), durations.end(), [deadline](double us) { return us > deadline; });
+  const double mean = std::accumulate(durations.begin(), durations.end(), 0.0) / durations.size();
+  std::cout << taps << ',' << callback_size << ',' << std::fixed << std::setprecision(2) << mean << ','
+            << percentile(durations, 0.50) << ',' << percentile(durations, 0.95) << ',' << percentile(durations, 0.99)
+            << ',' << durations.back() << ',' << deadline << ',' << overruns << ',' << sink + cache_sink * 0.0f << '\n';
+}
+} // namespace
+
+int main(int argc, char** argv)
+{
+  const int sample_rate = 48000;
+  const double seconds = argc > 3 ? std::atof(argv[3]) : 2.0;
+  const nam::LinearImplementation implementation =
+    argc > 4 && std::string(argv[4]) == "direct" ? nam::LinearImplementation::Direct : nam::LinearImplementation::FFT;
+  const bool cold_cache = argc > 5 && std::string(argv[5]) == "cold";
+  const bool verify = argc > 5 && std::string(argv[5]) == "verify";
+  const std::vector<int> taps = argc > 1 ? std::vector<int>{std::atoi(argv[1])}
+                                         : std::vector<int>{1024, 2048, 4096, 8192, 48000, 240000, 1200000, 2880000};
+  const std::vector<int> callbacks = argc > 2 ? std::vector<int>{std::atoi(argv[2])} : std::vector<int>{32, 64, 128};
+
+  if (verify)
+    return verify_impulse_response(taps.front(), callbacks.front(), sample_rate) ? 0 : 1;
+
+  std::cout << "taps,callback,mean_us,p50_us,p95_us,p99_us,max_us,deadline_us,overruns,sink\n";
+  for (const int tap_count : taps)
+    for (const int callback_size : callbacks)
+      run(tap_count, callback_size, sample_rate, seconds, implementation, cold_cache);
+  return 0;
+}

tools/run_tests.cpp

@@ -100,6 +100,8 @@ int main()
   test_linear::test_direct_known_values();
   test_linear::test_fft_matches_direct_irregular_chunks();
   test_linear::test_auto_selection();
+  test_linear::test_fft_dispatch_table();
+  test_linear::test_fft_impulse_response_across_dispatch_sizes();
   test_linear::test_parse_implementation();
   test_linear::test_direct_process_realtime_safe();
   test_linear::test_fft_process_realtime_safe();

tools/test/test_linear.cpp

@@ -140,17 +140,45 @@ void test_auto_selection()
   assert(short_model.GetRequestedImplementation() == nam::LinearImplementation::Auto);
   assert(short_model.GetActiveImplementation() == nam::LinearImplementation::Direct);
 
-  const auto cutoff_weights = make_weights(256, false);
-  nam::Linear cutoff_model(1, 1, 256, false, cutoff_weights, 48000.0);
+  const auto cutoff_weights = make_weights(1024, false);
+  nam::Linear cutoff_model(1, 1, 1024, false, cutoff_weights, 48000.0);
   assert(cutoff_model.GetRequestedImplementation() == nam::LinearImplementation::Auto);
   assert(cutoff_model.GetActiveImplementation() == nam::LinearImplementation::Direct);
 
-  const auto fft_weights = make_weights(512, false);
-  nam::Linear fft_model(1, 1, 512, false, fft_weights, 48000.0);
+  const auto fft_weights = make_weights(2048, false);
+  nam::Linear fft_model(1, 1, 2048, false, fft_weights, 48000.0);
   assert(fft_model.GetRequestedImplementation() == nam::LinearImplementation::Auto);
   assert(fft_model.GetActiveImplementation() == nam::LinearImplementation::FFT);
 }
 
+void test_fft_dispatch_table()
+{
+  assert(nam::linear::select_implementation(1024) == nam::LinearImplementation::Direct);
+  assert(nam::linear::select_implementation(1025) == nam::LinearImplementation::FFT);
+  assert(nam::linear::select_fft_plan(1024).direct_taps == 128);
+  assert(nam::linear::select_fft_plan(8192).max_partition_size == 2048);
+  assert(nam::linear::select_fft_plan(48000).max_partition_size == 4096);
+  assert(nam::linear::select_fft_plan(240000).max_partition_size == 8192);
+  assert(nam::linear::select_fft_plan(2880000).max_partition_size == 8192);
+}
+
+void test_fft_impulse_response_across_dispatch_sizes()
+{
+  const std::vector<int> receptive_fields{1024, 2048, 4096, 8192, 48000};
+  for (const int receptive_field : receptive_fields)
+  {
+    const auto weights = make_weights(receptive_field, false);
+    nam::Linear model(1, 1, receptive_field, false, weights, 48000.0, nam::LinearImplementation::FFT);
+    std::vector<NAM_SAMPLE> input(receptive_field + 257, 0.0);
+    input[0] = 1.0;
+    const auto output = process_model(model, input, {1, 17, 32, 63, 128, 511});
+    for (int i = 0; i < receptive_field; ++i)
+      assert_near(output[i], weights[i], 5.0e-5);
+    for (size_t i = receptive_field; i < output.size(); ++i)
+      assert_near(output[i], 0.0, 5.0e-5);
+  }
+}
+
 void test_parse_implementation()
 {
   assert(nam::linear::parse_implementation("auto") == nam::LinearImplementation::Auto);