score-avnd-amuencha/Amuencha/SpiralDisplay.cpp

#include "SpiralDisplay.hpp"

Amuencha::SpiralDisplay::SpiralDisplay()
    : min_midi_note{24}
    , max_midi_note{72}
{
  for (int i{0}; i < 12; i++) note_positions[i] = polar(.80f, half_pi - i * two_pi / 12);
}

void Amuencha::SpiralDisplay::compute_frequencies()
{
  // Now the spiral
  // 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 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);
  // one more bound than number of bins
  display_bins.resize(num_bins + 1);
  bin_sizes.resize(num_bins);
  spiral_positions.resize(num_bins + 1);
  spiral_r_a.resize(num_bins + 1);
  const float rmin = 0.1;
  const float rmax = 0.9;
  // The spiral and bounds are the same independently of how
  // the log space is divided into notes (e.g. 12ET)
  // Make it so c is on the y axis. Turn clockwise because people are
  // 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;

  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.);
    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] = polar(spiral_r_a[b].r, spiral_r_a[b].a);
  }

  // repeat one more time to avoid a second for loops
  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;
  spiral_r_a[num_bins].a = theta_min + (theta_max - theta_min) * bratio;
  spiral_positions[num_bins] = polar(spiral_r_a[num_bins].r, spiral_r_a[num_bins].a);

  for (int b{0}; b < num_bins; ++b)
    bin_sizes[b] = display_bins[b + 1] - display_bins[b];

  for (int id{0}; id < num_ID; ++id)
  {
    display_spectrum[id].resize(num_bins);
    fill(display_spectrum[id].begin(), display_spectrum[id].end(), 0.);
  }
}

void Amuencha::SpiralDisplay::power_handler(int ID, const std::vector<float> &reassigned_frequencies, const std::vector<float> &power_spectrum)
{
  fill(display_spectrum[ID].begin(), display_spectrum[ID].end(), 0.);

  // simple histogram-like sum, assuming power entries are normalized
  int nidx = reassigned_frequencies.size();
  for (int idx{0}; idx < nidx; ++idx)
  {
    float rf = reassigned_frequencies[idx];
    int ri = idx;
    // reassigned frequencies are never too far off the original
    //if (rf>display_bins[idx] && rf<display_bins[idx+1]) ri = idx;
    //else...
    while (rf<display_bins[ri])
    {
      --ri;
      if (ri==-1) break;
    }
    if (ri==-1) continue; // ignore this frequency, it is below display min
    while (rf>display_bins[ri+1])
    {
      ++ri;
      if (ri==nidx) break;
    }
    if (ri==nidx) continue; // ignore this frequency, it is above display max

    // Normalization:
    // - for a given frequency, the sine/window size dependency was already
    //   handled in the frequency analyzer
    // - but the result should not depend on how many frequencies are provided:
    //   increasing the resolution should not increase the power
    // => we need a kind of density, not just the histogram-like sum of powers
    //    falling into each bin
    // - consider the energy is coming from all the original bin size & sum
    // - This way, using finer bins do not increase the total sum
    display_spectrum[ID][ri] += power_spectrum[idx] * bin_sizes[idx];
  }

  // - Then, spread on the destination bin for getting uniform density
  //   measure independently of the target bin size
  for (int idx{0}; idx < nidx; ++idx) display_spectrum[ID][idx] /= bin_sizes[idx];
}