From 29283ede6340ca7525e47041b9689c3372d2a74c Mon Sep 17 00:00:00 2001 From: thibaudk Date: Tue, 29 Oct 2024 23:32:02 +0000 Subject: [PATCH] [Analyzer] use original frequency analyzer directly --- Amuencha/AmuenchaModel.cpp | 241 +------------------------------------ Amuencha/AmuenchaModel.hpp | 63 ++++------ Amuencha/AmuenchaUi.hpp | 7 +- Amuencha/SpiralDisplay.cpp | 2 +- Amuencha/SpiralDisplay.hpp | 2 - CMakeLists.txt | 5 +- 6 files changed, 32 insertions(+), 288 deletions(-) diff --git a/Amuencha/AmuenchaModel.cpp b/Amuencha/AmuenchaModel.cpp index 66c863d..df09fe4 100644 --- a/Amuencha/AmuenchaModel.cpp +++ b/Amuencha/AmuenchaModel.cpp @@ -1,249 +1,14 @@ #include "Amuencha.hpp" -#include "sse_mathfun.h" - -#include -#include - -#include -#include -#include namespace Amuencha { -void Analyser::prepare(halp::setup info) +void Model::prepare(halp::setup info) { sampling_rate = info.rate; - samplerate_div_2pi = sampling_rate/two_pi; } -void Analyser::operator()(halp::tick t) +void Model::operator()(halp::tick t) { - using namespace std; - - if (chunks.empty()) - chunks.emplace_back(make_pair(inputs.audio[0], t.frames)); - else - chunks[0] = make_pair(inputs.audio[0], t.frames); - - int new_data_size = 0; - for (auto c: chunks) new_data_size += c.second; - // shift the old data to make room for the new - int new_data_pos = big_buffer.size() - new_data_size; - //cout << "data " << new_data_size << " " << big_buffer.size() << endl; - // std::copy can cope with overlapping regions in this copy direction - if (new_data_pos > 0) - copy(big_buffer.begin() + new_data_size, - big_buffer.end(), - big_buffer.begin()); - // now copy each chunk at its position in the big buffer - for (auto c: chunks) - { - if (new_data_pos < 0) - { - // discard too old chunks - if (c.second <= -new_data_pos) - { - new_data_pos += c.second; - continue; - } - // partial copy of chunks that fall on the edge - copy(c.first + new_data_pos, c.first + c.second, &big_buffer[0]); - new_data_pos = c.second + new_data_pos; - continue; - } - - copy(c.first, c.first + c.second, &big_buffer[new_data_pos]); - new_data_pos += c.second; - } - - // Apply the filter bank - float *bbend = &big_buffer[0] + big_buffer.size(); - for (int idx{0}; idx < frequencies.size(); ++idx) - { - const auto& ws = windowed_sines[idx]; - int wsize = ws.size(); - float* sig = bbend - wsize; - v4sf acc = {0.f, 0.f, 0.f, 0.f}; - for (int i{0}; i < wsize; ++i) acc += ws[i] * sig[i]; - float norm = acc[0] * acc[0] + acc[1]*acc[1]; - float reassign = frequencies[idx]; - if (norm > 0) - reassign -= (acc[0] * acc[3] - acc[1] * acc[2]) * samplerate_div_2pi / norm; - - reassigned_frequencies[idx] = reassign; - power_spectrum[idx] = norm * power_normalization_factors[idx]; - } - - send_message({.min = inputs.min, - .max = inputs.max, - .reassigned_frequencies = this->reassigned_frequencies, - .power_spectrum = this->power_spectrum}); -} - -void Analyser::analyzer_setup(float max_buffer_duration) -{ - using namespace boost::math::float_constants; - using namespace std; - - // Prepare the windows - //std::vector> windowed_sines; - windowed_sines.resize(frequencies.size()); - power_normalization_factors.resize(frequencies.size()); - - int big_buffer_size = 0; - - for (int idx{0}; idx < frequencies.size(); ++idx) - { - // for each freq, span at least 20 periods for more precise measurements - // This still gives reasonable latencies, e.g. 50ms at 400Hz, 100ms at 200Hz, 400ms at 50Hz... - // Could also span more for even better measurements, with larger - // computation cost and latency - float f = frequencies[idx]; - int window_size = (int)(min(inputs.periods / f, max_buffer_duration * 0.001f) - * sampling_rate); - - vector window(window_size); - vector window_deriv(window_size); - - if (!read_from_cache(window, window_deriv)) - { - initialize_window(window); - initialize_window_deriv(window_deriv); - write_to_cache(window, window_deriv); - } - windowed_sines[idx].resize(window_size); - float wsum = 0; - - for (int i{0}; i < window_size;) - { - if (i < window_size - 4) - { - v4sf tfs = { - (float)(i - window_size - 1) / sampling_rate, - (float)(i + 1 - window_size - 1) / sampling_rate, - (float)(i + 2 - window_size - 1) / sampling_rate, - (float)(i + 3 - window_size - 1) / sampling_rate - }; - tfs *= (float)(-two_pi*f); - v4sf sin_tf, cos_tf; - sincos_ps(tfs, &sin_tf, &cos_tf); - - for (int j{0}; j < 3; ++j) - { - v4sf ws = { - cos_tf[j] * window[i + j], - sin_tf[j] * window[i + j], - cos_tf[j] * window_deriv[i + j], - sin_tf[j] * window_deriv[i + j] - }; - windowed_sines[idx][i + j] = ws; - wsum += window[i + j]; - } - i += 4; - continue; - } - float t = (float)(i - window_size - 1) / sampling_rate; - float re = cosf(-two_pi * t * f); - float im = sinf(-two_pi * t * f); - v4sf ws = { - re * window[i], - im * window[i], - re * window_deriv[i], - im * window_deriv[i] - }; - windowed_sines[idx][i] = ws; - wsum += window[i]; - ++i; - } - power_normalization_factors[idx] = 1. / (wsum * wsum); - big_buffer_size = max(big_buffer_size, window_size); - } - - big_buffer.clear(); - // fill with 0 signal content to start with - big_buffer.resize(big_buffer_size, 0.f); -} - -void Analyser::initialize_window(std::vector& window) -{ - // Kaiser window with a parameter of alpha=3 that nullifies the window on edges - int size = window.size(); - const float two_over_N = 2. / size; - const float alpha = 3.; - const float alpha_pi = alpha * pi; - const float inv_denom = 1. / boost::math::cyl_bessel_i(0., alpha_pi); - for (int i{0}; i < size; ++i) - { - float p = i * two_over_N - 1.; - window[i] = boost::math::cyl_bessel_i(0., alpha_pi * sqrt(1. - p * p)) * inv_denom; - } -} - -void Analyser::initialize_window_deriv(std::vector& window) -{ - // Derivative of the Kaiser window with a parameter of alpha=3 that nullifies the window on edges - int size = window.size(); - const float two_over_N = 2. / size; - const float alpha = 3.; - const float alpha_pi = alpha * pi; - const float inv_denom = 1. / boost::math::cyl_bessel_i(0., alpha_pi); - for (int i{1}; i < size; ++i) - { - float p = i * two_over_N - 1.; - window[i] = boost::math::cyl_bessel_i(1., alpha_pi * sqrt(1. - p * p)) * - inv_denom * alpha_pi / sqrt(1. - p * p) * (-p) * two_over_N; - } - // lim I1(x)/x as x->0 = 1/2 - window[0] = 0.5 * inv_denom * alpha_pi * alpha_pi * two_over_N; -} - -bool Analyser::read_from_cache(std::vector &window, std::vector &window_deriv) -{ - auto it = mem_win_cache.find(window.size()); - - if (it != mem_win_cache.end()) - { - window = it->second; - auto itd = mem_winderiv_cache.find(window.size()); - - if (itd != mem_winderiv_cache.end()) - { - window_deriv = itd->second; - return true; - } - // else, load from disk - } - - using namespace std; - // TODO: make the cache location parametrizable (and an option to not use it) - ifstream file(".amuencha_cache/w"+to_string(window.size())+".bin", ios::binary); - if (file.is_open()) - { - file.read(reinterpret_cast(&window[0]), window.size()*sizeof(float)); - file.read(reinterpret_cast(&window_deriv[0]), window_deriv.size()*sizeof(float)); - - if (file.tellg() != (window.size() + window_deriv.size()) * sizeof(float)) - { - cerr << "Error: invalid cache .amuencha_cache/w"+to_string(window.size())+".bin\n"; - return false; - } - return true; - } - return false; -} - -void Analyser::write_to_cache(std::vector &window, std::vector &window_deriv) -{ - using namespace std; -#if defined(_WIN32) || defined(_WIN64) - mkdir(".amuencha_cache"); -#else - mkdir(".amuencha_cache", 0755); -#endif - ofstream file(".amuencha_cache/w"+to_string(window.size())+".bin", ios::binary|ios::trunc); - file.write(reinterpret_cast(&window[0]), window.size()*sizeof(float)); - file.write(reinterpret_cast(&window_deriv[0]), window_deriv.size()*sizeof(float)); - mem_win_cache[window.size()] = window; - mem_winderiv_cache[window.size()] = window_deriv; + analyzer.new_data(inputs.audio[0], t.frames); } } diff --git a/Amuencha/AmuenchaModel.hpp b/Amuencha/AmuenchaModel.hpp index a4e006c..c2b9703 100644 --- a/Amuencha/AmuenchaModel.hpp +++ b/Amuencha/AmuenchaModel.hpp @@ -5,9 +5,11 @@ #include #include +#include "FrequencyAnalyzer.hpp" + namespace Amuencha { -class Analyser +class Model { public: halp_meta(name, "Amuencha") @@ -33,16 +35,16 @@ public: halp::fixed_audio_bus<"Input", float, 1> audio; struct : halp::hslider_i32<"Min", halp::range{.min = 0, .max = 127, .init = 24}> { - void update(Analyser& self) + void update(Model& self) { - self.send_message({.min = this->value, .max = self.inputs.max.value}); + self.send_message({.min = this->value, .max = self.inputs.max}); } } min; struct : halp::hslider_i32<"Min", halp::range{.min = 0, .max = 127, .init = 72}> { - void update(Analyser& self) + void update(Model& self) { - self.send_message({.min = self.inputs.min.value, .max = this->value}); + self.send_message({.min = self.inputs.min, .max = this->value}); } } max; halp::spinbox_i32<"Periods", @@ -51,10 +53,19 @@ public: void process_message(const std::vector& frequencies) { - this->frequencies = frequencies; - reassigned_frequencies = frequencies; - power_spectrum.resize(frequencies.size()); - analyzer_setup(); + analyzer.setup(sampling_rate, + frequencies, + [&] (const std::vector& r_f, + const std::vector& p_s) + { + send_message({.min = inputs.min, + .max = inputs.max, + .reassigned_frequencies = r_f, + .power_spectrum = p_s}); + }, + inputs.periods); + + analyzer.start(QThread::NormalPriority); } struct outs @@ -74,41 +85,9 @@ public: // UI is defined in another file to keep things clear. struct ui; - // max_buffer_duration in milliseconds, specifies the largest buffer for computing the frequency content - // At lower frequencies, long buffers are needed for accurate frequency separation. - // When that max buffer duration is reached, then it is capped and the frequency resolution decreases - // Too low buffers also limit the min_freq, duration must be >= period - void analyzer_setup(float max_buffer_duration = 500); - private: - // new data chunks arrived since the last periodic processing - std::vector> chunks; - - // The window is the usual Kaiser with alpha=3 - static void initialize_window(std::vector& window); - static void initialize_window_deriv(std::vector& window); - // The filter bank. One filter per frequency - // The 4 entries in the v4sf are the real, imaginary parts of the windowed - // sine wavelet, and the real, imaginary parts of the derived windowed sine - // used for reassigning the power spectrum. - // Hopefully, with SIMD, computing all 4 of them is the same price as just one - // TODO: v8sf and compute 2 freqs at the same time - typedef float v4sf __attribute__ ((vector_size (16))); - std::vector> windowed_sines; - std::vector frequencies; - std::vector power_normalization_factors; float sampling_rate; - float samplerate_div_2pi; - std::vector big_buffer; - std::vector reassigned_frequencies; - std::vector power_spectrum; - - // caching computations for faster init - // on disk for persistence between executions, - // in memory for avoiding reloading from disk when changing the spiral size - bool read_from_cache(std::vector& window, std::vector& window_deriv); - void write_to_cache(std::vector& window, std::vector& window_deriv); - std::map> mem_win_cache, mem_winderiv_cache; + FrequencyAnalyzer analyzer; }; } diff --git a/Amuencha/AmuenchaUi.hpp b/Amuencha/AmuenchaUi.hpp index 2a7228a..95b646a 100644 --- a/Amuencha/AmuenchaUi.hpp +++ b/Amuencha/AmuenchaUi.hpp @@ -1,12 +1,13 @@ #pragma once -#include #include + +#include #include "SpiralDisplay.hpp" namespace Amuencha { -struct Analyser::ui +struct Model::ui { using enum halp::colors; using enum halp::layouts; @@ -43,7 +44,7 @@ struct Analyser::ui // Receive a message on the UI thread from the processing thread static void process_message(ui& self, const processor_to_ui& msg) { - self.spiral.set_min_max_notes(msg.min, msg.max); + // self.spiral.set_min_max_notes(msg.min, msg.max); if (msg.power_spectrum.empty() || msg.reassigned_frequencies.empty()) diff --git a/Amuencha/SpiralDisplay.cpp b/Amuencha/SpiralDisplay.cpp index 3e7b1dc..8a5f15d 100644 --- a/Amuencha/SpiralDisplay.cpp +++ b/Amuencha/SpiralDisplay.cpp @@ -94,7 +94,7 @@ void Amuencha::SpiralDisplay::compute_frequencies() void Amuencha::SpiralDisplay::power_handler(const std::vector& reassigned_frequencies, const std::vector& power_spectrum) { - // fill(display_spectrum.begin(), display_spectrum.end(), 0.); + fill(display_spectrum.begin(), display_spectrum.end(), 0.); // simple histogram-like sum, assuming power entries are normalized int nidx = reassigned_frequencies.size(); diff --git a/Amuencha/SpiralDisplay.hpp b/Amuencha/SpiralDisplay.hpp index 3c9afd8..927502d 100644 --- a/Amuencha/SpiralDisplay.hpp +++ b/Amuencha/SpiralDisplay.hpp @@ -30,8 +30,6 @@ struct SpiralDisplay if (display_bins.empty()) compute_frequencies(); - ctx.set_stroke_color({0, 0, 0, 255}); - for (int i{0}; i < 12; i++) { ctx.move_to(half, half); diff --git a/CMakeLists.txt b/CMakeLists.txt index 9bf1086..ac32af9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,9 +23,10 @@ avnd_score_plugin_add( Amuencha/AmuenchaUi.hpp Amuencha/SpiralDisplay.hpp Amuencha/SpiralDisplay.cpp - Amuencha/sse_mathfun.h + Amuencha/FrequencyAnalyzer.hpp + Amuencha/FrequencyAnalyzer.cpp TARGET amuencha - MAIN_CLASS Analyser + MAIN_CLASS Model NAMESPACE Amuencha )