Guest User

Untitled

a guest
Jun 28th, 2020
90
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Dart 7.42 KB | None | 0 0
  1. // Public domain license
  2.   double coefficientOfDetermination(List<double> yObserved, List<double> yModel) {
  3.     // See https://en.wikipedia.org/wiki/Coefficient_of_determination
  4.     if (yModel.length != yObserved.length || yModel.length < 2) throw ArgumentError();
  5.     // Sum of squared residuals
  6.     double ssRes = zip([yObserved, yModel]).fold(0, (a, x) => a + math.pow(x[0] - x[1], 2));
  7.     double yObservedMean = yObserved.fold(0, (a, x) => a + x) / yObserved.length;
  8.     double ssTot = yObserved.fold(0, (a, x) => a + math.pow(x - yObservedMean, 2));
  9.     return 1 - (ssRes / ssTot);
  10.   }
  11.  
  12.   double SSRes(List<double> yObserved, List<double> yModel) {
  13.     if (yModel.length != yObserved.length || yModel.length < 2) throw ArgumentError();
  14.     // Sum of squared residuals
  15.     return zip([yObserved, yModel]).fold(0, (a, x) => a + math.pow(x[0] - x[1], 2));
  16.   }
  17.  
  18.   double sampleFreq(int k) { // IIRC the two tiers are to save CPU cycles
  19.     if (k <= 100) return 0.01 * k; // 0-1hz
  20.     else return 1 + 0.25 * (k - 100); // >1hz
  21.   }
  22.  
  23.   double getFreq(List<double> ts, List<double> xs) { // IIRC Non-uniform DFT
  24.       int N = xs.length;
  25.       List<double> X = List<double>.generate(N+100, (x) => 0.0);
  26.       List<double> XReal = List<double>.generate(N+100, (x) => 0.0);
  27.       List<double> XImag = List<double>.generate(N+100, (x) => 0.0);
  28.       double windowInterval = ts[ts.length - 1] - ts[0];
  29.       double freqUnit = (windowInterval/1000);
  30.       List<double> p = ts.map((it) => (it - ts[0]) / (windowInterval)).toList();
  31.       List<double> f = List<double>.generate(100, (int it) => sampleFreq(it) * freqUnit);
  32.       double max = 0;
  33.       double maxFreq = 0;
  34.       for (var k = 0; k < f.length; k += 1) {
  35.         for (var n = 0; n < N; n += 1) {
  36.           var realPart = xs[n] * math.cos(2 * math.pi * p[n] * f[k]);
  37.           var iPart = xs[n] * -1 * math.sin(2 * math.pi * p[n] * f[k]);
  38.           XReal[k] += realPart;
  39.           XImag[k] += iPart;
  40.         }
  41.         X[k] = math.sqrt(math.pow(XReal[k], 2) + math.pow(XImag[k], 2));
  42.         if (X[k] > max) {
  43.           max = X[k];
  44.           maxFreq = f[k]*2*math.pi/(windowInterval);
  45.       //maxPhase = atan(-XImag[k]/XReal[k])
  46.         }
  47.       }
  48.       return maxFreq;
  49.   }
  50.  
  51.   List<double> dft(List<double> ts, List<double> xs) { // IIRC this is a rudimentary vestigiality
  52.       int N = xs.length;
  53.       List<double> X = List<double>.generate(N+100, (x) => 0.0);
  54.       List<double> XReal = List<double>.generate(N+100, (x) => 0.0);
  55.       List<double> XImag = List<double>.generate(N+100, (x) => 0.0);
  56.       double windowInterval = ts[ts.length - 1] - ts[0];
  57.       double freqUnit = (windowInterval/1000);
  58.       List<double> p = ts.map((it) => (it - ts[0]) / (windowInterval)).toList();
  59.       List<double> f = List<double>.generate(100, (int it) => sampleFreq(it) * freqUnit);
  60.       for (var k = 0; k < f.length; k += 1) {
  61.         for (var n = 0; n < N; n += 1) {
  62.           var realPart = 0.01 * xs[n] * math.cos(2 * math.pi * p[n] * f[k]);
  63.           var iPart = 0.01 * xs[n] * -1 * math.sin(2 * math.pi * p[n] * f[k]);
  64.           XReal[k] += realPart;
  65.           XImag[k] += iPart;
  66.         }
  67.         X[k] = math.sqrt(math.pow(XReal[k], 2) + math.pow(XImag[k], 2));
  68.       }
  69.       return X;
  70.   }
  71.  
  72.   List<double> sine(List<double> t, double m, double a, double p, double f) {
  73.     return t.map((t) => m + a*math.cos(p+t*f)).toList();
  74.   }
  75.  
  76.   double gridSearchPhase(List<double> dataWindow, List<double> timeWindow, double m) {
  77.     double rMin = double.maxFinite;
  78.     double pMin = 0;
  79.  
  80.     for (int i = 0; i < 50; i += 1) {
  81.       var param = 2*3.141596 * (i/50);
  82.       List<double> yModel = sine(timeWindow, m, _amp, param, _freq);
  83.       double r = SSRes(dataWindow, yModel);
  84.       if (r < rMin) {
  85.         rMin = r;
  86.         //debugPrint(rMax.toString());
  87.         pMin = param;
  88.       }
  89.     }
  90.     return pMin;
  91.     //return new math.Random().nextDouble() * 6.28;
  92.   }
  93.  
  94.   double gridSearchAmp(List<double> dataWindow, List<double> timeWindow, double m) {
  95.     double ampMax = (dataWindow.reduce(math.max) - dataWindow.reduce(math.min))/2;
  96.     double rMin = double.maxFinite;
  97.     double pMin = 0;
  98.  
  99.     for (int i = 0; i < 50; i += 1) {
  100.       var param = ampMax * 0.5 + ampMax * 0.5 * (i/50);
  101.       List<double> yModel = sine(timeWindow, m, param, _phase, _freq);
  102.       double r = SSRes(dataWindow, yModel);
  103.       if (r < rMin) {
  104.         rMin = r;
  105.         //debugPrint(rMax.toString());
  106.         pMin = param;
  107.       }
  108.     }
  109.     return pMin;
  110.     //return new math.Random().nextDouble() * 6.28;
  111.   }
  112.  
  113.  
  114.   List<double> _fitSineToData(List<double> dataWindow, List<double> timeWindow) {
  115.     double sampleMean = dataWindow.reduce((a, x) => a + x) / dataWindow.length;
  116.     //double threshold = median(List<double>.from(dataWindow.map((x) => (x - sampleMedian).abs())));
  117.     List<double> dataSansBaseline = dataWindow.map((x) => x - sampleMean).toList();
  118.     double estFreq = getFreq(timeWindow, dataSansBaseline);
  119.     //platform.invokeMethod('fitSineToData', <String, dynamic>{
  120.     //  //'xs': Float64List.fromList(List<double>.generate(180, (i) => i.toDouble()-1)),
  121.     // 'xs': Float64List.fromList(timeWindow),
  122.     //  //'ys': Float64List.fromList(dataWindow.map((x) => math.max((x - sampleMedian).abs(), threshold)*(x-sampleMedian).sign).toList())
  123.     // 'ys': Float64List.fromList(dataWindow.map((x) => (x - sampleMean)).toList())
  124.     //}).then(
  125.     //  (result) => setEstimates(result),
  126.     //  onError: (e) => debugPrint("error")
  127.     //);
  128.     _amp = ((dataWindow.reduce(math.max) - dataWindow.reduce(math.min))/2)*0.9;
  129.     _freq = estFreq;
  130.     _phase = gridSearchPhase(dataWindow, timeWindow, sampleMean);
  131.     _amp = gridSearchAmp(dataWindow, timeWindow, sampleMean);
  132.     var yModel = sine(timeWindow, sampleMean, _amp, _phase, _freq);
  133.     //var yModel = List<double>.generate(dataWindow.length, (i) => sampleMean + _amp * math.cos((_phase)+timeWindow[i]*_freq));
  134.     _goodnessOfFit = coefficientOfDetermination(dataWindow, yModel);
  135.     return yModel;
  136.   }
  137.  
  138.   void setEstimates(Float64List estimates) {
  139.     _amp = estimates[0];
  140.     _freq = estimates[1];
  141.     _phase = estimates[2];
  142.   }
  143.  
  144.   double _goodnessOfFit = 0.0;
  145.   double _amp = 0.0;
  146.   double _freq = 0.0;
  147.   double _phase = 0.0;
  148.  
  149.  
  150. This is free and unencumbered software released into the public domain.
  151.  
  152. Anyone is free to copy, modify, publish, use, compile, sell, or
  153. distribute this software, either in source code form or as a compiled
  154. binary, for any purpose, commercial or non-commercial, and by any
  155. means.
  156.  
  157. In jurisdictions that recognize copyright laws, the author or authors
  158. of this software dedicate any and all copyright interest in the
  159. software to the public domain. We make this dedication for the benefit
  160. of the public at large and to the detriment of our heirs and
  161. successors. We intend this dedication to be an overt act of
  162. relinquishment in perpetuity of all present and future rights to this
  163. software under copyright law.
  164.  
  165. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  166. EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  167. MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  168. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  169. OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  170. ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  171. OTHER DEALINGS IN THE SOFTWARE.
  172.  
  173. For more information, please refer to <http://unlicense.org/>
Advertisement
Add Comment
Please, Sign In to add comment