Fixed bug in implementation.

Author: arkottke

source/WangRathjePeakCalculator.cpp

@@ -42,34 +42,54 @@ double WangRathjePeakCalculator::calcDurationRms(
                 duration, oscFreq, oscDamping, siteTransFunc);
 
     bool validTransFunc = false;
-    validTransFunc &= siteTransFunc.size() == _freqs.size();
-    // Make sure there are some values greater than 1
-    if (validTransFunc) {
+    // Make sure there are some values greater than 1.1
+    if (siteTransFunc.size() == _freqs.size()) {
         for (std::complex<double> stf : siteTransFunc) {
-            if (abs(abs(stf) - 1.) > 1E-1) {
+            if (abs(stf) > 1.1) {
                 validTransFunc = true;
                 break;
             }
         }
     }
+    qDebug() << validTransFunc;
 
     // Modify the duration for the soil response
     if (validTransFunc) {
+        // Compute the expected rock oscillator duration
+
+        // Equation 4a
+        const double fLim = 5.274 * std::pow(duration, -0.640);
+        double ratio = 1;
+        // Equation 2
+        if (0.1 <= oscFreq && oscFreq < fLim) {
+            // Equation 4b
+            const double dur0 = 31.858 * std::pow(duration, -0.849);
+            // Equation 4c
+            const double durMin = 1.009 * duration / (3.583 + duration);
+            // Equation 3b
+            const double b = 1 / (dur0 - durMin);
+            // Equation 3a
+            const double a = (1 / (dur0 - 1) - b) * (fLim - 0.1);
+            // Equation 2
+            ratio = (dur0 - (oscFreq - 0.1) / (a + b * (oscFreq - 0.1)));
+        }
+
+        const double durOscRock = ratio * duration;
+
+        // Compute the expected soil oscillator duration
         QList<double> modeFreqs;
         QList<double> modeAmps;
 
-        QVector<double> absSiteTransFunc;
-        for (std::complex<double> stf : siteTransFunc) {
-            absSiteTransFunc << abs(stf);
+        QVector<double> absSiteTransFunc(siteTransFunc.size());
+        for (int i = 0; i < absSiteTransFunc.size(); ++i) {
+            absSiteTransFunc[i] = std::abs(siteTransFunc.at(i));
         }
 
         // Compute the location of the peaks in the transfer function
         const int offset = 2;
         bool valid;
         for (int i = offset; i < (_freqs.size() - offset); ++i) {
             valid = true;
-            qDebug() << _freqs.at(i) << absSiteTransFunc.at(i);
-
             for (int j = (i - offset); j <= (i + offset); ++j) {
                 if (j < i) {
                     valid &= absSiteTransFunc.at(j) < absSiteTransFunc.at(i);
@@ -90,19 +110,21 @@ double WangRathjePeakCalculator::calcDurationRms(
         // Amplitude / frequency ratio of the first mode
         double afRatio = modeAmps.at(0) / modeFreqs.at(0);
 
-        double a;
+        double incrMax;
         double c;
         double m;
         double incr = 0;
         for (int i = 0; i < modeFreqs.size(); ++i) {
             c = _coefs.at(i).a * afRatio + _coefs.at(i).b * std::pow(afRatio, 2);
             m = _coefs.at(i).d * afRatio + _coefs.at(i).e * std::pow(afRatio, 2);
-            a = c * std::exp(-duration / m);
-            incr += (a * std::exp(-std::pow(std::log(oscFreq / modeFreqs.at(i)), 2) /
+            incrMax = c * std::exp(-duration / m);
+            incr += (incrMax * std::exp(-std::pow(std::log(oscFreq / modeFreqs.at(i)), 2) /
                     (2 * std::pow(_coefs.at(i).sd, 2)))
             );
         }
-        durationRms += incr;
+        const double durOscSoil = durOscRock + incr;
+        qDebug() << oscFreq << durationRms << durationRms * (durOscSoil/ durOscRock);
+        durationRms *= (durOscSoil / durOscRock);
     }
     return durationRms;
 }