Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- // Public domain license
- double coefficientOfDetermination(List<double> yObserved, List<double> yModel) {
- // See https://en.wikipedia.org/wiki/Coefficient_of_determination
- if (yModel.length != yObserved.length || yModel.length < 2) throw ArgumentError();
- // Sum of squared residuals
- double ssRes = zip([yObserved, yModel]).fold(0, (a, x) => a + math.pow(x[0] - x[1], 2));
- double yObservedMean = yObserved.fold(0, (a, x) => a + x) / yObserved.length;
- double ssTot = yObserved.fold(0, (a, x) => a + math.pow(x - yObservedMean, 2));
- return 1 - (ssRes / ssTot);
- }
- double SSRes(List<double> yObserved, List<double> yModel) {
- if (yModel.length != yObserved.length || yModel.length < 2) throw ArgumentError();
- // Sum of squared residuals
- return zip([yObserved, yModel]).fold(0, (a, x) => a + math.pow(x[0] - x[1], 2));
- }
- double sampleFreq(int k) { // IIRC the two tiers are to save CPU cycles
- if (k <= 100) return 0.01 * k; // 0-1hz
- else return 1 + 0.25 * (k - 100); // >1hz
- }
- double getFreq(List<double> ts, List<double> xs) { // IIRC Non-uniform DFT
- int N = xs.length;
- List<double> X = List<double>.generate(N+100, (x) => 0.0);
- List<double> XReal = List<double>.generate(N+100, (x) => 0.0);
- List<double> XImag = List<double>.generate(N+100, (x) => 0.0);
- double windowInterval = ts[ts.length - 1] - ts[0];
- double freqUnit = (windowInterval/1000);
- List<double> p = ts.map((it) => (it - ts[0]) / (windowInterval)).toList();
- List<double> f = List<double>.generate(100, (int it) => sampleFreq(it) * freqUnit);
- double max = 0;
- double maxFreq = 0;
- for (var k = 0; k < f.length; k += 1) {
- for (var n = 0; n < N; n += 1) {
- var realPart = xs[n] * math.cos(2 * math.pi * p[n] * f[k]);
- var iPart = xs[n] * -1 * math.sin(2 * math.pi * p[n] * f[k]);
- XReal[k] += realPart;
- XImag[k] += iPart;
- }
- X[k] = math.sqrt(math.pow(XReal[k], 2) + math.pow(XImag[k], 2));
- if (X[k] > max) {
- max = X[k];
- maxFreq = f[k]*2*math.pi/(windowInterval);
- //maxPhase = atan(-XImag[k]/XReal[k])
- }
- }
- return maxFreq;
- }
- List<double> dft(List<double> ts, List<double> xs) { // IIRC this is a rudimentary vestigiality
- int N = xs.length;
- List<double> X = List<double>.generate(N+100, (x) => 0.0);
- List<double> XReal = List<double>.generate(N+100, (x) => 0.0);
- List<double> XImag = List<double>.generate(N+100, (x) => 0.0);
- double windowInterval = ts[ts.length - 1] - ts[0];
- double freqUnit = (windowInterval/1000);
- List<double> p = ts.map((it) => (it - ts[0]) / (windowInterval)).toList();
- List<double> f = List<double>.generate(100, (int it) => sampleFreq(it) * freqUnit);
- for (var k = 0; k < f.length; k += 1) {
- for (var n = 0; n < N; n += 1) {
- var realPart = 0.01 * xs[n] * math.cos(2 * math.pi * p[n] * f[k]);
- var iPart = 0.01 * xs[n] * -1 * math.sin(2 * math.pi * p[n] * f[k]);
- XReal[k] += realPart;
- XImag[k] += iPart;
- }
- X[k] = math.sqrt(math.pow(XReal[k], 2) + math.pow(XImag[k], 2));
- }
- return X;
- }
- List<double> sine(List<double> t, double m, double a, double p, double f) {
- return t.map((t) => m + a*math.cos(p+t*f)).toList();
- }
- double gridSearchPhase(List<double> dataWindow, List<double> timeWindow, double m) {
- double rMin = double.maxFinite;
- double pMin = 0;
- for (int i = 0; i < 50; i += 1) {
- var param = 2*3.141596 * (i/50);
- List<double> yModel = sine(timeWindow, m, _amp, param, _freq);
- double r = SSRes(dataWindow, yModel);
- if (r < rMin) {
- rMin = r;
- //debugPrint(rMax.toString());
- pMin = param;
- }
- }
- return pMin;
- //return new math.Random().nextDouble() * 6.28;
- }
- double gridSearchAmp(List<double> dataWindow, List<double> timeWindow, double m) {
- double ampMax = (dataWindow.reduce(math.max) - dataWindow.reduce(math.min))/2;
- double rMin = double.maxFinite;
- double pMin = 0;
- for (int i = 0; i < 50; i += 1) {
- var param = ampMax * 0.5 + ampMax * 0.5 * (i/50);
- List<double> yModel = sine(timeWindow, m, param, _phase, _freq);
- double r = SSRes(dataWindow, yModel);
- if (r < rMin) {
- rMin = r;
- //debugPrint(rMax.toString());
- pMin = param;
- }
- }
- return pMin;
- //return new math.Random().nextDouble() * 6.28;
- }
- List<double> _fitSineToData(List<double> dataWindow, List<double> timeWindow) {
- double sampleMean = dataWindow.reduce((a, x) => a + x) / dataWindow.length;
- //double threshold = median(List<double>.from(dataWindow.map((x) => (x - sampleMedian).abs())));
- List<double> dataSansBaseline = dataWindow.map((x) => x - sampleMean).toList();
- double estFreq = getFreq(timeWindow, dataSansBaseline);
- //platform.invokeMethod('fitSineToData', <String, dynamic>{
- // //'xs': Float64List.fromList(List<double>.generate(180, (i) => i.toDouble()-1)),
- // 'xs': Float64List.fromList(timeWindow),
- // //'ys': Float64List.fromList(dataWindow.map((x) => math.max((x - sampleMedian).abs(), threshold)*(x-sampleMedian).sign).toList())
- // 'ys': Float64List.fromList(dataWindow.map((x) => (x - sampleMean)).toList())
- //}).then(
- // (result) => setEstimates(result),
- // onError: (e) => debugPrint("error")
- //);
- _amp = ((dataWindow.reduce(math.max) - dataWindow.reduce(math.min))/2)*0.9;
- _freq = estFreq;
- _phase = gridSearchPhase(dataWindow, timeWindow, sampleMean);
- _amp = gridSearchAmp(dataWindow, timeWindow, sampleMean);
- var yModel = sine(timeWindow, sampleMean, _amp, _phase, _freq);
- //var yModel = List<double>.generate(dataWindow.length, (i) => sampleMean + _amp * math.cos((_phase)+timeWindow[i]*_freq));
- _goodnessOfFit = coefficientOfDetermination(dataWindow, yModel);
- return yModel;
- }
- void setEstimates(Float64List estimates) {
- _amp = estimates[0];
- _freq = estimates[1];
- _phase = estimates[2];
- }
- double _goodnessOfFit = 0.0;
- double _amp = 0.0;
- double _freq = 0.0;
- double _phase = 0.0;
- This is free and unencumbered software released into the public domain.
- Anyone is free to copy, modify, publish, use, compile, sell, or
- distribute this software, either in source code form or as a compiled
- binary, for any purpose, commercial or non-commercial, and by any
- means.
- In jurisdictions that recognize copyright laws, the author or authors
- of this software dedicate any and all copyright interest in the
- software to the public domain. We make this dedication for the benefit
- of the public at large and to the detriment of our heirs and
- successors. We intend this dedication to be an overt act of
- relinquishment in perpetuity of all present and future rights to this
- software under copyright law.
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
- EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
- OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
- ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
- OTHER DEALINGS IN THE SOFTWARE.
- For more information, please refer to <http://unlicense.org/>
Advertisement
Add Comment
Please, Sign In to add comment