+// it's perhaps not ideal to _first_ find the peak and _then_ the harmonics --
+// ideally, one should find the set of all peaks and then determine the likely
+// base from that... something for later. :-)
+std::pair<double, double> adjust_for_overtones(std::pair<double, double> base, double *in, unsigned num_samples)
+{
+ double mu = base.first, var = 1.0 / (base.second * base.second);
+
+ for (unsigned i = 2; i < 10; ++i) {
+ unsigned middle = unsigned(floor(freq_to_bin(base.first, num_samples) * i + 0.5));
+ unsigned lower = middle - (i+1)/2, upper = middle + (i+1)/2;
+
+ if (upper >= num_samples)
+ upper = num_samples - 2;
+
+ // printf("Searching in [%u,%u] = %f..%f\n", lower, upper, bin_to_freq(lower, num_samples), bin_to_freq(upper, num_samples));
+
+ // search for a peak in this interval
+ double best_harmonics_freq = -1.0;
+ double best_harmonics_amp = -1.0;
+ for (unsigned j = lower; j <= upper; ++j) {
+ if (in[j] > in[j-1] && in[j] > in[j+1]) {
+ std::pair<double, double> peak =
+ interpolate_peak(in[j - 1],
+ in[j],
+ in[j + 1]);
+
+ if (peak.second > best_harmonics_amp) {
+ best_harmonics_freq = bin_to_freq(j + peak.first, num_samples);
+ best_harmonics_amp = peak.second;
+ }
+ }
+ }
+
+ if (best_harmonics_amp <= 0.0)
+ continue;
+
+ //printf("Found overtone %u at %.2f (amp=%5.2fdB)\n", i, best_harmonics_freq,
+ // best_harmonics_amp);
+
+ double this_mu = best_harmonics_freq / double(i);
+ double this_var = 1.0 / (best_harmonics_amp * best_harmonics_amp);
+
+ double k = var / (var + this_var);
+ mu = (1.0 - k) * mu + k * this_mu;
+ var *= (1.0 - k);
+ }
+ return std::make_pair(mu, base.second);
+}
+