thibaud.keller/score-avnd-amuencha
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Compare commits

merge into: thibaud.keller:master
Branches Tags
thibaud.keller:master
thibaud.keller:Worker
thibaud.keller:QThread
...
pull from: thibaud.keller:QThread
Branches Tags
thibaud.keller:Worker
thibaud.keller:QThread
thibaud.keller:master
Sign in to create a new pull request.

5 commits
master ... QThread

Author SHA1 Message Date
thibaudk
cba01d2519 [Readme] update TODO list 2024-12-06 11:00:14 +00:00
thibaudk
c8e051a4a3 [frequencies] seperate compute_frequencies betwin analyzer and display 2024-11-01 15:35:00 +00:00
thibaudk
a4a4ba249e [ui] add CustomSlider 2024-11-01 00:01:20 +00:00
thibaudk
eaf8314a68 [display] first working extremly dirty poc 2024-10-30 13:07:46 +00:00
thibaudk
29283ede63 [Analyzer] use original frequency analyzer directly 2024-10-29 23:32:02 +00:00
13 changed files with 630 additions and 356 deletions
Show stats Expand all files Collapse all files

2
Amuencha/Amuencha.hpp
View file

@ -1,4 +1,4 @@
#pragma once
#include <Amuencha/AmuenchaModel.hpp>
#include <Amuencha/AmuenchaUi.hpp>
#include <Amuencha/AmuenchaUi.hpp>

253
Amuencha/AmuenchaModel.cpp
View file

@ -1,249 +1,24 @@
#include "Amuencha.hpp"
#include "sse_mathfun.h"
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include <sys/stat.h>
#include <iostream>
#include <fstream>
namespace Amuencha
{
void Analyser::prepare(halp::setup info)
void Model::prepare(setup info)
{
sampling_rate = info.rate;
samplerate_div_2pi = sampling_rate/two_pi;
if (analyzer.isRunning()) return;
analyzer.setup(info.rate,
[&] (const std::vector<float>& r_f,
const std::vector<float>& p_s)
{
this->send_message({.reassigned_frequencies = r_f,
.power_spectrum = p_s});
});
analyzer.start(QThread::NormalPriority);
}
void Analyser::operator()(halp::tick t)
void Model::operator()(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<std::vector<v4sf>> 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<float> window(window_size);
vector<float> 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<float>& 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<float>& 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<float> &window, std::vector<float> &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<char*>(&window[0]), window.size()*sizeof(float));
file.read(reinterpret_cast<char*>(&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<float> &window, std::vector<float> &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<char*>(&window[0]), window.size()*sizeof(float));
file.write(reinterpret_cast<char*>(&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);
}
}

74
Amuencha/AmuenchaModel.hpp
View file

@ -5,9 +5,11 @@
#include <halp/controls.hpp>
#include <halp/meta.hpp>
#include "FrequencyAnalyzer.hpp"
namespace Amuencha
{
class Analyser
class Model
{
public:
halp_meta(name, "Amuencha")
@ -15,11 +17,8 @@ public:
halp_meta(c_name, "amuencha")
halp_meta(uuid, "b37351b4-7b8d-4150-9e1c-708eec9182b2")
// This one will be memcpy'd as it is a trivial type
struct processor_to_ui
{
int min;
int max;
std::vector<float> reassigned_frequencies;
std::vector<float> power_spectrum;
};
@ -31,84 +30,29 @@ public:
struct ins
{
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)
{
self.send_message({.min = this->value, .max = self.inputs.max.value});
}
} min;
struct : halp::hslider_i32<"Min", halp::range{.min = 0, .max = 127, .init = 72}>
{
void update(Analyser& self)
{
self.send_message({.min = self.inputs.min.value, .max = this->value});
}
} max;
halp::spinbox_i32<"Periods",
halp::range{.min = 0, .max = 99, .init = 30}> periods;
halp::hslider_i32<"Min", halp::range{.min = 0, .max = 127, .init = 24}> min;
halp::hslider_i32<"Max", halp::range{.min = 0, .max = 127, .init = 72}> max;
halp::spinbox_i32<"Periods", halp::range{.min = 0, .max = 99, .init = 30}> periods;
} inputs;
void process_message(const std::vector<float>& frequencies)
{
this->frequencies = frequencies;
reassigned_frequencies = frequencies;
power_spectrum.resize(frequencies.size());
analyzer_setup();
}
struct outs
{
halp::midi_bus<"Output"> midi;
} outputs;
using setup = halp::setup;
void prepare(halp::setup info);
void prepare(setup info);
// Do our processing for N samples
using tick = halp::tick;
// Defined in the .cpp
void operator()(halp::tick t);
void operator()(tick t);
// 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<std::pair<float*,int>> chunks;
// The window is the usual Kaiser with alpha=3
static void initialize_window(std::vector<float>& window);
static void initialize_window_deriv(std::vector<float>& 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<std::vector<v4sf>> windowed_sines;
std::vector<float> frequencies;
std::vector<float> power_normalization_factors;
float sampling_rate;
float samplerate_div_2pi;
std::vector<float> big_buffer;
std::vector<float> reassigned_frequencies;
std::vector<float> 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<float>& window, std::vector<float>& window_deriv);
void write_to_cache(std::vector<float>& window, std::vector<float>& window_deriv);
std::map<int, std::vector<float>> mem_win_cache, mem_winderiv_cache;
FrequencyAnalyzer analyzer;
};
}

23
Amuencha/AmuenchaUi.hpp
View file

@ -1,12 +1,14 @@
#pragma once
#include <Amuencha/AmuenchaModel.hpp>
#include <halp/layout.hpp>
#include <Amuencha/AmuenchaModel.hpp>
#include "SpiralDisplay.hpp"
#include "CustomSlider.hpp"
namespace Amuencha
{
struct Analyser::ui
struct Model::ui
{
using enum halp::colors;
using enum halp::layouts;
@ -18,33 +20,32 @@ struct Analyser::ui
struct
{
halp_meta(layout, vbox)
halp::item<&ins::min> min;
halp::item<&ins::max> max;
halp::custom_control<CustomSlider, &ins::min> min;
halp::custom_control<CustomSlider, &ins::max> max;
halp::item<&ins::periods> periods;
} controls;
halp::custom_actions_item<SpiralDisplay> spiral{.x = 0, .y = 0};
// Define the communication between UI and processor.
struct bus
{
std::function<void(const std::vector<float>&)> send_message;
// Set up connections
void init(ui& self)
{
self.spiral.on_new_frequencies = [&]
self.controls.min.on_changed = [&] (int min)
{
send_message(self.spiral.get_frequencies());
self.spiral.set_min_max_notes(min, self.controls.max.value);
};
self.controls.max.on_changed = [&] (int max)
{
self.spiral.set_min_max_notes(self.controls.min.value, max);
};
}
// 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);
if (msg.power_spectrum.empty() ||
msg.reassigned_frequencies.empty())
return;

2
Amuencha/CustomSlider.cpp Normal file
View file

@ -0,0 +1,2 @@
#include "CustomSlider.hpp"

74
Amuencha/CustomSlider.hpp Normal file
View file

@ -0,0 +1,74 @@
#pragma once
#include <halp/custom_widgets.hpp>
#include <avnd/wrappers/controls.hpp>
#include <avnd/concepts/painter.hpp>
struct CustomSlider
{
// Same as above
static constexpr double width() { return 100.; }
static constexpr double height() { return 20.; }
// Needed for changing the ui. It's the type above - it's already defined as-is
// in the helpers library.
halp::transaction<int> transaction;
// Called when the value changes from the host software.
void set_value(const auto& control, int value)
{
this->value = avnd::map_control_to_01(control, value / 127.f);
}
// When transaction.update() is called, this converts the value in the slider
// into one fit for the control definition passed as argument.
static auto value_to_control(auto& control, int value)
{
return avnd::map_control_from_01(control, value / 127.f);
}
// Paint method: same as above
void paint(avnd::painter auto ctx)
{
ctx.set_stroke_color({200, 200, 200, 255});
ctx.set_stroke_width(2.);
ctx.set_fill_color({120, 120, 120, 255});
ctx.begin_path();
ctx.draw_rect(0., 0., width(), height());
ctx.fill();
ctx.stroke();
ctx.begin_path();
ctx.set_fill_color({90, 90, 90, 255});
ctx.draw_rect(2., 2., (width() - 4) * (value / 127.f), (height() - 4));
ctx.fill();
}
// Return true to handle the event. x, y, are the positions of the item in local coordinates.
bool mouse_press(double x, double y)
{
transaction.start();
mouse_move(x, y);
return true;
}
// Obvious :-)
void mouse_move(double x, double y)
{
const int res{static_cast<int>(std::clamp(x / width(), 0., 1.) * 127)};
transaction.update(res);
on_changed(value);
}
// Same
void mouse_release(double x, double y)
{
mouse_move(x, y);
transaction.commit();
}
int value{};
std::function<void(int)> on_changed{[](int){}};
};

367
Amuencha/FrequencyAnalyzer.cpp Normal file
View file

@ -0,0 +1,367 @@
/*
Analyseur de MUsique et ENtraînement au CHAnt
This file is released under either of the two licenses below, your choice:
- LGPL v2.1 or later, https://www.gnu.org
The GNU Lesser General Public Licence, version 2.1 or,
at your option, any later version.
- CeCILL-C, http://www.cecill.info
The CeCILL-C license is more adapted to the French laws,
but can be converted to the GNU LGPL.
You can use, modify and/or redistribute the software under the terms of any
of these licences, which should have been provided to you together with this
sofware. If that is not the case, you can find a copy of the licences on
the indicated web sites.
By Nicolas . Brodu @ Inria . fr
See http://nicolas.brodu.net/programmation/amuencha/ for more information
*/
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/constants/constants.hpp>
#include <sys/stat.h>
#include "sse_mathfun.h"
#include "FrequencyAnalyzer.hpp"
using namespace std;
using namespace boost::math::float_constants;
Amuencha::FrequencyAnalyzer::FrequencyAnalyzer(QObject *parent) : QThread(parent)
{
status = NO_DATA;
}
Amuencha::FrequencyAnalyzer::~FrequencyAnalyzer()
{
mutex.lock();
status = QUIT_NOW;
// Do not even wait the end of a cycle, quit
condition.wakeOne();
mutex.unlock();
wait(CYCLE_PERIOD*2);
terminate();
wait(); // until run terminates
}
void Amuencha::FrequencyAnalyzer::new_data(float* chunk, int size)
{
// producer, called from another thread.
// Must not wake the main thread uselessly to notify that new data is available
// The idea is that the new_data may be called very fast from a RT audio thread
// but that frequencies are computed at most 10 times/second (or whatever the cyclic period is set)
// => decouple thread frequencies
// And when no data is here (recording off, no song...), the main thread is fully passive (no CPU hog)
mutex.lock();
// Store a pointer here and not a full data copy, which is much faster
// The real data is stored in the AudioRecording object
// Kind of duplicates the AudioRecording list, however:
// - This way, there is no need for mutex/lock in the AudioRecording structure
// - The position of the last unprocessed chunk within that structure need not be stored there
chunks.emplace_back(make_pair(chunk, size));
// set the flag that data has arrived
status = HAS_NEW_DATA;
// IF AND ONLY IF the thread was blocked forever, then wake it up
if (waiting_time != CYCLE_PERIOD) {
// and now resume the cyclic scheduling
waiting_time = CYCLE_PERIOD;
condition.wakeOne();
}
// Otherwise, do NOT wake the other thread, keep the low-freq cycle to decrease load
mutex.unlock();
}
void Amuencha::FrequencyAnalyzer::run()
{
// Solution with a wait condition + time
// other possible solution = timer, but that would need to be stopped
mutex.lock();
// waiting_time is a mutex-protected info
waiting_time = CYCLE_PERIOD;
// loop starts with mutex locked
while (true) {
condition.wait(&mutex, waiting_time);
if (status==QUIT_NOW) break;
if (status==HAS_NEW_DATA) {
// Swap the chunks to a local variable, so:
// - the class chunks becomes empty
// - the audio thread can feed it more data while computing frequencies
vector<pair<float*,int>> chunks;
chunks.swap(this->chunks);
status = NO_DATA; // will be updated if new data indeed arrives
mutex.unlock();
// Now, we can take the time to do the frequency computations
data_mutex.lock();
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];
}
// Notify our listener that new power/frequency content has arrived
power_handler(reassigned_frequencies, power_spectrum);
// setup can now lock and change data structures if needed
data_mutex.unlock();
// relock for the condition wait
mutex.lock();
continue;
}
// No more data ? Force waiting until data arrives
waiting_time = ULONG_MAX;
// keep the lock for next loop
}
mutex.unlock();
}
void Amuencha::FrequencyAnalyzer::setup(float sampling_rate,
PowerHandler&& handler,
int min_midi_note,
int max_midi_note,
float periods,
float max_buffer_duration)
{
// Block data processing while changing the data structures
data_mutex.lock();
this->samplerate_div_2pi = sampling_rate/two_pi;
// Start with A440, but this could be parametrizable as well
const float fref = 440;
const float log2_fref = log2(fref);
const int aref = 69; // use the midi numbering scheme, because why not
float log2_fmin = (min_midi_note - aref) / 12. + log2_fref;
float log2_fmax = (max_midi_note - aref) / 12. + log2_fref;
int num_bins = max_midi_note - min_midi_note;
frequencies.resize(num_bins);
for (int b{0}; b < num_bins; ++b)
{
float bratio = (float)b / (num_bins - 1.);
frequencies[b] = exp2(log2_fmin + (log2_fmax - log2_fmin) * bratio);
}
this->reassigned_frequencies = frequencies;
this->power_spectrum.resize(frequencies.size());
power_handler = std::move(handler);
// Prepare the windows
//std::vector<std::vector<v4sf>> 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(periods / f, max_buffer_duration * 0.001f) * sampling_rate);
vector<float> window(window_size);
vector<float> 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);
// Processing can resume with the new data structures in place.
data_mutex.unlock();
}
void Amuencha::FrequencyAnalyzer::invalidate_samples()
{
mutex.lock();
data_mutex.lock();
chunks.clear();
data_mutex.unlock();
mutex.unlock();
}
void Amuencha::FrequencyAnalyzer::initialize_window(std::vector<float>& 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 Amuencha::FrequencyAnalyzer::initialize_window_deriv(std::vector<float>& 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 Amuencha::FrequencyAnalyzer::read_from_cache(std::vector<float>& window, std::vector<float>& 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
}
// 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<char*>(&window[0]), window.size() * sizeof(float));
file.read(reinterpret_cast<char*>(&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 Amuencha::FrequencyAnalyzer::write_to_cache(std::vector<float>& window, std::vector<float>& window_deriv)
{
#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<char*>(&window[0]), window.size() * sizeof(float));
file.write(reinterpret_cast<char*>(&window_deriv[0]), window_deriv.size() * sizeof(float));
mem_win_cache[window.size()] = window;
mem_winderiv_cache[window.size()] = window_deriv;
}

117
Amuencha/FrequencyAnalyzer.hpp Normal file
View file

@ -0,0 +1,117 @@
/*
Analyseur de MUsique et ENtraînement au CHAnt
This file is released under either of the two licenses below, your choice:
- LGPL v2.1 or later, https://www.gnu.org
The GNU Lesser General Public Licence, version 2.1 or,
at your option, any later version.
- CeCILL-C, http://www.cecill.info
The CeCILL-C license is more adapted to the French laws,
but can be converted to the GNU LGPL.
You can use, modify and/or redistribute the software under the terms of any
of these licences, which should have been provided to you together with this
sofware. If that is not the case, you can find a copy of the licences on
the indicated web sites.
By Nicolas . Brodu @ Inria . fr
See http://nicolas.brodu.net/programmation/amuencha/ for more information
*/
#ifndef AUDIOINPUTTHREAD_H
#define AUDIOINPUTTHREAD_H
#include <QThread>
#include <QMutex>
#include <QWaitCondition>
#include <vector>
#include <map>
#include <functional>
namespace Amuencha
{
class FrequencyAnalyzer : public QThread
{
public:
FrequencyAnalyzer(QObject *parent = 0);
~FrequencyAnalyzer();
// called by the RT audio thread to feed new data
void new_data(float *chunk, int size);
// Arguments are: frequency bins [f,f+1), and power in each bin
// hence the first vector size is 1 more than the second
// TODO if useful: make a proper listener API with id, etc
typedef std::function<void(const std::vector<float>&,const std::vector<float>&)> PowerHandler;
PowerHandler power_handler;
// sampling rate in Hz
// freqs in Hz
// PowerHandler for the callback
// 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 setup(float sampling_rate,
PowerHandler&& handler,
int min_midi_note = 24,
int max_midi_note = 72,
float periods = 20,
float max_buffer_duration = 500);
// call to remove all existing chunk references
// this may cause signal loss, but this is usually called precisely when the signal is lost...
void invalidate_samples();
protected:
void run() override;
// Multi-threading related variables
static const unsigned long CYCLE_PERIOD = 20; // in milliseconds
QMutex mutex, data_mutex;
QWaitCondition condition;
unsigned long waiting_time = CYCLE_PERIOD;
enum Status : std::uint8_t
{
NO_DATA = 0,
HAS_NEW_DATA = 1,
QUIT_NOW = 2
};
Status status;
// new data chunks arrived since the last periodic processing
std::vector<std::pair<float*,int>> chunks;
// The window is the usual Kaiser with alpha=3
static void initialize_window(std::vector<float>& window);
static void initialize_window_deriv(std::vector<float>& 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<std::vector<v4sf>> windowed_sines;
std::vector<float> frequencies;
std::vector<float> power_normalization_factors;
float samplerate_div_2pi;
std::vector<float> big_buffer;
std::vector<float> reassigned_frequencies;
std::vector<float> 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<float>& window, std::vector<float>& window_deriv);
void write_to_cache(std::vector<float>& window, std::vector<float>& window_deriv);
std::map<int, std::vector<float>> mem_win_cache, mem_winderiv_cache;
};
} // namespace Amuencha
#endif // AUDIOINPUTTHREAD_H

28
Amuencha/SpiralDisplay.cpp
View file

@ -5,7 +5,6 @@ Amuencha::SpiralDisplay::SpiralDisplay()
, max_midi_note{72}
, gain{1.f}
, visual_fading{1}
, on_new_frequencies{[]{}}
{
for (int i{0}; i < 12; i++)
note_positions[i] = std::polar(.9f, half_pi - i * two_pi / 12);
@ -24,11 +23,6 @@ void Amuencha::SpiralDisplay::set_min_max_notes(int min_midi_note, int max_midi_
// update();
}
std::vector<float> Amuencha::SpiralDisplay::get_frequencies() const noexcept
{
return frequencies;
}
void Amuencha::SpiralDisplay::compute_frequencies()
{
// Now the spiral
@ -38,16 +32,17 @@ void Amuencha::SpiralDisplay::compute_frequencies()
const int aref = 69; // use the midi numbering scheme, because why not
float log2_fmin = (min_midi_note - aref) / 12. + log2_fref;
float log2_fmax = (max_midi_note - aref) / 12. + log2_fref;
int approx_pix_bin_width = 3;
// int approx_pix_bin_width = 3;
// number of frequency bins is the number of pixels
// along the spiral path / approx_pix_bin_width
// According to mathworld, the correct formula for the path length
// from the origin involves sqrt and log computations.
// Here, we just want some approximate pixel count
// => use all circles for the approx
int num_octaves = (max_midi_note - min_midi_note + 11) / 12;
float approx_num_pix = 0.5 * half * pi * num_octaves;
int num_bins = (int)(approx_num_pix / approx_pix_bin_width);
// int num_octaves = (max_midi_note - min_midi_note + 11) / 12;
// float approx_num_pix = 0.5 * half * pi * num_octaves;
// int num_bins = (int)(approx_num_pix / approx_pix_bin_width);
int num_bins = max_midi_note - min_midi_note;
// one more bound than number of bins
display_bins.resize(num_bins + 1);
bin_sizes.resize(num_bins);
@ -61,21 +56,18 @@ void Amuencha::SpiralDisplay::compute_frequencies()
// used to it (e.g. wikipedia note circle)
const float theta_min = half_pi - two_pi * (min_midi_note % 12) / 12;
// wrap in anti-trigonometric direction
const float theta_max = theta_min - two_pi * (max_midi_note - min_midi_note) / 12;
const float theta_max = theta_min - two_pi * (max_midi_note - min_midi_note - 1) / 12;
frequencies.resize(num_bins);
for (int b{0}; b < num_bins; ++b)
{
float bratio = (float)b / (num_bins - 1.);
frequencies[b] = exp2(log2_fmin + (log2_fmax - log2_fmin) * bratio);
bratio = (float)(b - 0.5) / (float)(num_bins - 1.);
float bratio = (float)(b - 0.5) / (float)(num_bins - 1.);
display_bins[b] = exp2(log2_fmin + (log2_fmax - log2_fmin) * bratio);
spiral_r_a[b].r = rmin + (rmax - rmin) * bratio;
spiral_r_a[b].a = theta_min + (theta_max - theta_min) * bratio;
spiral_positions[b] = std::polar(spiral_r_a[b].r, spiral_r_a[b].a);
}
// repeat one more time to avoid a second for loops
// repeat one more time to avoid a second for loop
float bratio = (float)(num_bins - 0.5) / (float)(num_bins - 1.);
display_bins[num_bins] = exp2(log2_fmin + (log2_fmax - log2_fmin) * bratio);
spiral_r_a[num_bins].r = rmin + (rmax - rmin) * bratio;
@ -87,14 +79,12 @@ void Amuencha::SpiralDisplay::compute_frequencies()
display_spectrum.resize(num_bins);
fill(display_spectrum.begin(), display_spectrum.end(), 0.);
on_new_frequencies();
}
void Amuencha::SpiralDisplay::power_handler(const std::vector<float>& reassigned_frequencies,
const std::vector<float>& 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();

22
Amuencha/SpiralDisplay.hpp
View file

@ -4,7 +4,7 @@
#include <avnd/concepts/painter.hpp>
#include <array>
#include <ranges>
#include <vector>
#include <complex>
namespace Amuencha
@ -30,8 +30,8 @@ struct SpiralDisplay
if (display_bins.empty())
compute_frequencies();
ctx.set_stroke_color({0, 0, 0, 255});
// Draw note axis and names
// TODO : replace with CSV background
for (int i{0}; i < 12; i++)
{
ctx.move_to(half, half);
@ -60,7 +60,7 @@ struct SpiralDisplay
for (int b{0}; b < display_spectrum.size(); ++b)
{
float amplitude = 0.8 / num_octaves * std::min(1.f, display_spectrum[b] * gain);
float amplitude = 0.8 / num_octaves * std::min(1.f, display_spectrum[b] * 100);
//if (display_spectrum[b]>0) cout << display_spectrum[b] << endl;
// power normalised between 0 and 1 => 0.1 = spiral branch
float r = spiral_r_a[b].r + amplitude;
@ -71,22 +71,19 @@ struct SpiralDisplay
ctx.line_to(x(p.real()), y(p.imag()));
}
for (int b = spiral_positions.size() - 1; b >= 0; --b)
ctx.line_to(x(spiral_positions[b].real()), y(spiral_positions[b].imag()));
// for (int b = spiral_positions.size() - 1; b >= 0; --b)
// ctx.line_to(x(spiral_positions[b].real()), y(spiral_positions[b].imag()));
ctx.stroke();
ctx.update();
}
[[nodiscard]] std::vector<float> get_frequencies() const noexcept;
std::function<void()> on_new_frequencies;
// // Callback when the power spectrum is available at the prescribed frequencies
// // The ID is that of the caller, setting the color of the display
void power_handler(const std::vector<float>& reassigned_frequencies,
const std::vector<float>& power_spectrum);
void compute_frequencies();
private:
static const constexpr std::string_view note_names[12]
{"C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B"};
@ -106,9 +103,6 @@ private:
int min_midi_note, max_midi_note, visual_fading;
// central frequencies (log space)
std::vector<float> frequencies;
// local copy for maintaining the display, adapted to the drawing bins
std::vector<float> display_spectrum;
@ -135,8 +129,6 @@ private:
{
return half - y * half;
}
void compute_frequencies();
};
}

7
CMakeLists.txt
View file

@ -20,12 +20,15 @@ avnd_score_plugin_add(
Amuencha/Amuencha.hpp
Amuencha/AmuenchaModel.hpp
Amuencha/AmuenchaModel.cpp
Amuencha/FrequencyAnalyzer.hpp
Amuencha/FrequencyAnalyzer.cpp
Amuencha/AmuenchaUi.hpp
Amuencha/SpiralDisplay.hpp
Amuencha/SpiralDisplay.cpp
Amuencha/sse_mathfun.h
Amuencha/CustomSlider.hpp
Amuencha/CustomSlider.cpp
TARGET amuencha
MAIN_CLASS Analyser
MAIN_CLASS Model
NAMESPACE Amuencha
)

17
README.md
View file

@ -2,16 +2,25 @@
An [ossia score](https://ossia.io) add-on cloned from the [avendish template](https://github.com/ossia-templates/score-avnd-simple-template)\
It aims at porting Nicolas Brodu's [frequency analyzer](https://nicolas.brodu.net/en/programmation/amuencha/index.html) to Avendish
![](amuencha.png)
## Features
* Draw background spiral and notes axis
## TODO
* Visualisation
- [ ] Redraw with min and max midi value
- [ ] Rimplement all draw function without storing painters
* Ui
- [X] Redraw with Min and Max midi value
- [X] Rimplement draw function without storing painters
- [ ] add an SVG backgound with fixed Axis and not names
- [ ] Min and Max note in one ranger rather than 2 sliders
- [ ] Display note Names and frequencies
* Port Backend analyzer
- [ ] Reimplement with avendish [threaded worker](https://celtera.github.io/avendish/advanced/workers.html)
- [ ] Asses the need for setting params therough [shared instance](https://gothub.ducks.party/celtera/avendish/blob/6fb720259c176442db32aafad67e27b5fe4b7b39/include/halp/shared_instance.hpp)
* Additional features
- [ ] Midi output
- [ ] audio output per frequency band
- [ ] MPE output
- [ ] Calculate pitch corelation between input channels (maybe draw from [S8-PDP](https://gothub.ducks.party/scrime-u-bordeaux/S8-PDP))

BIN
amuencha.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 65 KiB