Update WaveSpectrumEstimator.h

Author: mgrouch

src/spectrum/WaveSpectrumEstimator.h

@@ -132,96 +132,6 @@ class EIGEN_ALIGN_MAX WaveSpectrumEstimator {
         return (fp > 0.0) ? (1.0 / fp) : 0.0;
     }
 
-    PMFitResult fitPiersonMoskowitz() const {
-        Vec S_obs = lastSpectrum_;
-        for (int i = 0; i < Nfreq; ++i) {
-            if (S_obs[i] <= 0.0) S_obs[i] = 1e-12;
-        }
-
-        auto cost_fn = [&](double a, double fp) {
-            const double omega_p = 2.0 * M_PI * fp;
-            double cost = 0.0;
-            constexpr double beta = 0.74;
-
-            for (int i = 0; i < Nfreq - 1; ++i) {
-                const double df = df_[i];
-                const double f = freqs_[i];
-                double model = a * g * g * std::pow(2.0 * M_PI * f, -5.0)
-                             * std::exp(-beta * std::pow(omega_p / (2.0 * M_PI * f), 4.0));
-                if (model <= 0.0) model = 1e-12;
-                const double d = WaveSpectrumShared::safe_log(S_obs[i]) -
-                                 WaveSpectrumShared::safe_log(model);
-                cost += df * d * d;
-            }
-            return cost;
-        };
-
-        constexpr int N_fp_search = 32;
-        constexpr double fp_min = 0.05;
-        constexpr double fp_transition = 0.1;
-        constexpr double fp_max = 1.0;
-
-        std::array<double, N_fp_search> fp_grid{};
-        const int n_log = static_cast<int>(N_fp_search * 0.4);
-        const int n_lin = N_fp_search - n_log;
-
-        for (int i = 0; i < n_log; ++i) {
-            const double t = double(i) / double(n_log - 1);
-            fp_grid[i] = fp_min * std::pow(fp_transition / fp_min, t);
-        }
-        for (int i = 0; i < n_lin; ++i) {
-            const double t = double(i) / double(n_lin - 1);
-            fp_grid[n_log + i] = fp_transition + t * (fp_max - fp_transition);
-        }
-
-        double bestA = 1e-5;
-        double bestFp = fp_grid[0];
-        double bestC = std::numeric_limits<double>::infinity();
-
-        for (int ia = 0; ia < 8; ++ia) {
-            const double a = 1e-5 + ia * (1.0 - 1e-5) / 7.0;
-            for (int ifp = 0; ifp < N_fp_search; ++ifp) {
-                const double fp = fp_grid[ifp];
-                const double c = cost_fn(a, fp);
-                if (c < bestC) {
-                    bestC = c;
-                    bestA = a;
-                    bestFp = fp;
-                }
-            }
-        }
-
-        double alpha = bestA;
-        double fp = bestFp;
-        double stepA = 0.1;
-        double stepFp = 0.1;
-
-        for (int iter = 0; iter < 40; ++iter) {
-            bool improved = false;
-            double c = 0.0;
-
-            c = cost_fn(alpha + stepA, fp);
-            if (c < bestC) { bestC = c; alpha += stepA; improved = true; }
-
-            c = cost_fn(alpha - stepA, fp);
-            if (c < bestC) { bestC = c; alpha -= stepA; improved = true; }
-
-            c = cost_fn(alpha, fp + stepFp);
-            if (c < bestC) { bestC = c; fp += stepFp; improved = true; }
-
-            c = cost_fn(alpha, fp - stepFp);
-            if (c < bestC) { bestC = c; fp -= stepFp; improved = true; }
-
-            if (!improved) {
-                stepA *= 0.5;
-                stepFp *= 0.5;
-                if (stepA < 1e-12 && stepFp < 1e-12) break;
-            }
-        }
-
-        return {alpha, fp, bestC};
-    }
-
     void set_regularization_f0(double f0_hz) {
         reg_f0_hz = std::max(1e-6, f0_hz);
     }