~8 % improvement

not quite sure why this improved so much

changes are too complicated to trace it exactly

much of the improvement was related how I refactored the log10
calculations to be simpler and just log based

|Score type |MSE               |Min score         |Max score         |Mean score        |
|-----------|------------------|------------------|------------------|------------------|
|Zimtohrli  |0.090239608814795 |0.534520229950808 |0.772224479921073 |0.707988288758192 |
|ViSQOL     |0.115330916105424 |0.520833375452983 |0.801480831107469 |0.675101633981268 |
|2f         |0.129541391104905 |0.484687555319526 |0.797475783883375 |0.661870345773127 |
|PESQ       |0.147425552045669 |0.342342966279351 |0.841271127756762 |0.647128996775172 |
|CDPAM      |0.153471222942756 |0.441558428344727 |0.728779141125759 |0.620699318941738 |
|PARLAQ     |0.185057687192323 |0.445261140223642 |0.784370761057963 |0.587162756572532 |
|AQUA       |0.223207996944378 |0.331645933512413 |0.739286336419790 |0.547804951221731 |
|PEAQB      |0.225217321572038 |0.278744167467764 |0.851011116004117 |0.553935720513487 |
|DPAM       |0.315810440183130 |0.186717781679534 |0.690564701717118 |0.460415212267967 |
|WARP-Q     |0.339686211572685 |0.067600137543649 |0.777119464646524 |0.475793617709890 |
|GVPMOS     |0.412937133868407 |0.006851162794410 |0.783946603687895 |0.412912222208318 |

real	4m18.248s
user	221m18.093s
sys	46m19.418s

cpp/zimt/fourier_bank.cc

@@ -19,9 +19,87 @@
 namespace tabuli {
 
 float SimpleDb(float energy) {
-  static const float full_scale_sine_db = 78.26561963526045; // 78.3 ideally
-  static const float epsilon = 1.0033294789821357e-09;
-  return 10 * log10(energy + epsilon) + full_scale_sine_db;
+  // ideally 78.3 db
+  static const float full_scale_sine_db = 75.27901963526045;
+  static const float exp_full_scale_sine_db = exp(full_scale_sine_db);
+  // epsilon, but the biggest one you saw (~4.95e23)
+  static const float epsilon = 1.0033294789821357e-09 * exp_full_scale_sine_db;
+  // kMul allows faster log instead of log10 below, incorporating multiplying by 10 for decibel.
+  constexpr float kMul = 10.0/log(10);
+  return kMul * log(energy + epsilon);
+}
+
+void FinalizeDb(hwy::AlignedNDArray<float, 2>& channels, float mul,
+                size_t out_ix) {
+  double masker = 0.0;
+  static const double octaves_in_20_to_20000 = log(20000/20.)/log(2);
+  static const double octaves_per_rot =
+      octaves_in_20_to_20000 / float(kNumRotators - 1);
+  static const double masker_step_per_octave_up_0 = 15.892019717473835;
+  static const double masker_step_per_octave_up_1 = 21.852019717473834;
+  static const double masker_step_per_octave_up_2 = 20.79201971747383;
+  static const double masker_step_per_rot_up_0 = octaves_per_rot * masker_step_per_octave_up_0;
+  static const double masker_step_per_rot_up_1 = octaves_per_rot * masker_step_per_octave_up_1;
+  static const double masker_step_per_rot_up_2 = octaves_per_rot * masker_step_per_octave_up_2;
+  static const double masker_gap_up = 19.140338374861235;
+  static const float maskingStrengthUp = 0.1252262923615547;
+  static const float up_blur = 0.8738593591692092;
+  static const float fraction_up = 1.02;
+
+  static const double masker_step_per_octave_down = 42.33972783112732;
+  static const double masker_step_per_rot_down = octaves_per_rot * masker_step_per_octave_down;
+  static const double masker_gap_down = 19.66099875393617;
+  static const float maskingStrengthDown = 0.19329999999999992;
+  static const float down_blur = 0.714425315233319;
+
+  static const float min_limit = -11.397341001787765;
+  static const float fraction_down = 1.02;
+  // Scan frequencies from bottom to top, let lower frequencies to mask higher frequencies.
+  // 'masker' maintains the masking envelope from one bin to next.
+  for (int k = 0; k < kNumRotators; ++k) {
+    float v = SimpleDb(mul * channels[{out_ix}][k]);
+    if (v < min_limit) {
+      v = min_limit;
+    }
+    float v2 = (1 - up_blur) * v2 + up_blur * v;
+    if (k == 0) {
+      v2 = v;
+    }
+    if (masker < v2) {
+      masker = v2;
+    }
+    float mask = fraction_up * masker - masker_gap_up;
+    if (v < mask) {
+      v = maskingStrengthUp * mask + (1.0 - maskingStrengthUp) * v;
+    }
+    channels[{out_ix}][k] = v;
+    if (3 * k < kNumRotators) {
+      masker -= masker_step_per_rot_up_0;
+    } else if (3 * k < 2 * kNumRotators) {
+      masker -= masker_step_per_rot_up_1;
+    } else {
+      masker -= masker_step_per_rot_up_2;
+    }
+  }
+  // Scan frequencies from top to bottom, let higher frequencies to mask lower frequencies.
+  // 'masker' maintains the masking envelope from one bin to next.
+  masker = 0.0;
+  for (int k = kNumRotators - 1; k >= 0; --k) {
+    float v = channels[{out_ix}][k];
+    float v2 = (1 - down_blur) * v2 + down_blur * v;
+    if (k == kNumRotators - 1) {
+      v2 = v;
+    }
+    if (masker < v) {
+      masker = v;
+    }
+    float mask = fraction_down * masker - masker_gap_down;
+    if (v < mask) {
+      v = maskingStrengthDown * mask + (1.0 - maskingStrengthDown) * v;
+    }
+    channels[{out_ix}][k] = v;
+    masker -= masker_step_per_rot_down;
+  }
 }
 
 void Rotators::FilterAndDownsample(hwy::Span<const float> signal,
@@ -35,10 +113,7 @@ void Rotators::FilterAndDownsample(hwy::Span<const float> signal,
       int64_t input_ix = ii + zz;
       if (input_ix >= signal.size()) {
         if (out_ix < channels.shape()[0]) {
-          for (int k = 0; k < kNumRotators; ++k) {
-            channels[{out_ix}][k] =
-                SimpleDb(scaling_for_downsampling * channels[{out_ix}][k]);
-          }
+          FinalizeDb(channels, scaling_for_downsampling, out_ix);
         }
         if (out_ix != channels.shape()[0] - 1) {
           fprintf(stderr,
@@ -64,10 +139,7 @@ void Rotators::FilterAndDownsample(hwy::Span<const float> signal,
         }
       }
     }
-    for (int k = 0; k < kNumRotators; ++k) {
-      channels[{out_ix}][k] =
-          SimpleDb(scaling_for_downsampling * channels[{out_ix}][k]);
-    }
+    FinalizeDb(channels, scaling_for_downsampling, out_ix);
     ++out_ix;
     if (out_ix >= channels.shape()[0]) {
       return;
@@ -85,7 +157,7 @@ Rotators::Rotators(int num_channels, std::vector<float> frequency,
                    std::vector<float> filter_gains, const float sample_rate) {
   channel.resize(num_channels);
   static const float kWindow = 0.9996028710680265;
-  static const double kBandwidthMagic = 0.7322402492401391;
+  static const double kBandwidthMagic = 0.7328516996032982;
   for (int i = 0; i < kNumRotators; ++i) {
     // The parameter relates to the frequency shape overlap and window length
     // of triple leaking integrator.
@@ -95,7 +167,8 @@ Rotators::Rotators(int num_channels, std::vector<float> frequency,
     window[i] = std::pow(kWindow, bw * kBandwidthMagic);
     float windowM1 = 1.0f - window[i];
     float f = frequency[i] * 2.0f * M_PI / sample_rate;
-    const float gainer = 2.0f;
+    static const float full_scale_sine_db = exp(75.27901963526045);
+    const float gainer = 2.0f * sqrt(full_scale_sine_db);
     gain[i] = gainer * filter_gains[i] * pow(windowM1, 3.0);
     rot[0][i] = float(std::cos(f));
     rot[1][i] = float(-std::sin(f));

cpp/zimt/zimtohrli.cc

@@ -294,13 +294,8 @@ void Zimtohrli::Spectrogram(
   tabuli::Rotators rots(1, freqs, gains, cam_filterbank->sample_rate);
   rots.FilterAndDownsample(signal, energy_channels_db, downsample);
 
-  if (apply_masking) {
-    masking.CutFullyMasked(energy_channels_db, cam_filterbank->cam_delta,
-                           partial_energy_channels_db);
-  } else {
-    hwy::CopyBytes(energy_channels_db.data(), partial_energy_channels_db.data(),
-                   energy_channels_db.memory_size() * sizeof(float));
-  }
+  hwy::CopyBytes(energy_channels_db.data(), partial_energy_channels_db.data(),
+                 energy_channels_db.memory_size() * sizeof(float));
   if (apply_loudness) {
     loudness.PhonsFromSPL(partial_energy_channels_db,
                           cam_filterbank->thresholds_hz, spectrogram);