1 /* 2 * guiding_stats.h 3 * PHD2 Guiding 4 * 5 * Created by Bruce Waddington 6 * Copyright (c) 2018 Bruce Waddington 7 * All rights reserved. 8 * 9 * This source code is distributed under the following "BSD" license 10 * Redistribution and use in source and binary forms, with or without 11 * modification, are permitted provided that the following conditions are met: 12 * Redistributions of source code must retain the above copyright notice, 13 * this list of conditions and the following disclaimer. 14 * Redistributions in binary form must reproduce the above copyright notice, 15 * this list of conditions and the following disclaimer in the 16 * documentation and/or other materials provided with the distribution. 17 * Neither the name of Bret McKee, Dad Dog Development, 18 * Craig Stark, Stark Labs nor the names of its 19 * contributors may be used to endorse or promote products derived from 20 * this software without specific prior written permission. 21 * 22 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 23 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 25 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 26 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 27 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 28 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 29 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 30 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 31 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 32 * POSSIBILITY OF SUCH DAMAGE. 33 * 34 */ 35 36 #ifndef _GUIDING_STATS_H 37 #define _GUIDING_STATS_H 38 #include <deque> 39 40 // DescriptiveStats is used for basic statistics. Max, min, sigma and variance are computed on-the-fly as values are added to a dataset 41 // Applicable to any double values, no semantic assumptions made. Does not retain a list of values 42 class DescriptiveStats 43 { 44 private: 45 int count; 46 double runningS; // S denotes square of deltas from mean (variance) 47 double newS; 48 double runningMean; // current mean of dataset 49 double newMean; 50 double minValue; // current min value of dataset 51 double maxValue; // current max value of dataset 52 double lastValue; // for clients to easily compute deltas 53 double maxDelta; // max absolute (delta(value)) 54 55 public: 56 DescriptiveStats(); 57 ~DescriptiveStats(); 58 void AddValue(double Val); // Add a double value to the dataset 59 void ClearAll(); // Start over, reset all variables 60 unsigned int GetCount(); // Returns the count of the dataset 61 double GetLastValue(); // Returns the immediately previous value added to the dataset 62 double GetMean(); // Returns the mean value of the dataset 63 double GetSum(); // Sum of all added vars 64 double GetMinimum(); // Returns the min value 65 double GetMaximum(); // Returns the max value 66 double GetVariance(); // Variance for those who need it 67 double GetSigma(); // Returns the standard deviation 68 double GetPopulationSigma(); // Population sigma ('n' vs 'n-1') 69 double GetMaxDelta(); // Returns max of absolute delta(new - previous) values 70 }; 71 72 // High and Low pass filters can be used to filter a stream of data elements that can then be added to DescriptiveStats or AxisStats 73 // A low-pass filter will attenuate (dampen) high-frequency elements, a high-pass filter will do the opposite 74 // Examples: use a low-pass filter to emphasize low-frequency data fluctuations such as a slow linear drift 75 // use a high-pass filter to ignore linear drift but emphasize more rapid fluctuations 76 // Neither class retains any history of the data, the client is responsble for using the HPF/LPF values as needed 77 class HighPassFilter 78 { 79 private: 80 double alphaCutoff = 1.0; 81 double count; 82 double prevVal; 83 double hpfResult; 84 85 public: HighPassFilter()86 HighPassFilter() {}; 87 HighPassFilter(double CutoffPeriod, double SamplePeriod); 88 double AddValue(double NewVal); 89 double GetCurrentHPF(); 90 void Reset(); 91 }; 92 93 class LowPassFilter 94 { 95 private: 96 double alphaCutoff = 1.0; 97 double count; 98 double lpfResult; 99 100 public: LowPassFilter()101 LowPassFilter() {}; 102 LowPassFilter(double CutoffPeriod, double SamplePeriod); 103 double AddValue(double NewVal); 104 double GetCurrentLPF(); 105 void Reset(); 106 }; 107 108 // Support structure for use with AxisStats to keep a queue of guide star displacements and relative time values 109 // Timestamps are intended to be incremental, i.e seconds since start of guiding, and are used only for linear fit operations 110 struct StarDisplacement 111 { 112 double DeltaTime; 113 double StarPos; 114 bool Guided; 115 bool Reversal; 116 117 StarDisplacement(double When, double Where); 118 }; 119 120 // AxisStats and the StarDisplacement class can be used to collect and evaluate typical guiding data. Datasets can be windowed or not. 121 // Windowing means the data collection is limited to the most recent <n> entries. 122 // Windowed datasets will be automatically trimmed if AutoWindowSize > 0 or can be manually trimmed by client using RemoveOldestEntry() 123 class AxisStats 124 { 125 protected: 126 std::deque <StarDisplacement> guidingEntries; // queue of elements in dataset 127 unsigned int axisMoves; // number of times in window when guide pulse was non-zero 128 unsigned int axisReversals; // number of times in window when guide pulse caused a direction reversal 129 double prevMove; // value of guide pulse in next-to-last entry 130 double prevPosition; // value of guide star location in next-to-last entry 131 // Variables used to compute stats in windowed AxisStats 132 double sumX; // Sum of the x values (deltaT values) 133 double sumY; // Sum of the y values (star position) 134 double sumXY; // Sum of (x * y) 135 double sumXSq; // Sum of (x squared) 136 double sumYSq; // Sum of (y squared) 137 // Variables needed for windowed or non-windowed versions 138 double maxDisplacement; // maximum star position value in current dataset 139 double minDisplacement; // minimum star position value in current dataset 140 double maxDelta; // maximum absolute delta of incremental star deltas 141 int maxDeltaInx; 142 void InitializeScalars(); 143 144 public: 145 // Constructor for 3 types of instance: non-windowed, windowed with automatic trimming of size, windowed but with client controlling actual window size 146 AxisStats(); 147 ~AxisStats(); 148 149 // Add a guiding info element of relative time, guide star position, guide pulse amount 150 void AddGuideInfo(double DeltaT, double StarPos, double GuideAmt); 151 152 // Return a particular element from the current dataset 153 StarDisplacement GetEntry(unsigned int index) const; 154 155 void ClearAll(); 156 157 // Return the count of elements in the dataset 158 unsigned int GetCount() const; 159 160 // Get the last entry added to the dataset - useful if client needs to do difference operations for time, star position, or guide amount 161 StarDisplacement GetLastEntry() const; 162 163 // Get the maximum y value in the dataset 164 double GetMaxDisplacement() const; 165 166 // Get the minimum y value in the dataset 167 double GetMinDisplacement() const; 168 169 // Return stats for current dataset - min, max, sum, mean, variance, standard deviation (sigma) 170 double GetSum() const; 171 double GetMean() const; 172 double GetVariance() const; 173 double GetSigma() const; 174 double GetPopulationSigma() const; 175 double GetMedian() const; 176 double GetMaxDelta() const; 177 // Count of moves or reversals in current dataset 178 unsigned int GetMoveCount() const; 179 unsigned int GetReversalCount() const; 180 181 // Perform a linear fit on the star position values in the dataset. Return usual slope and y-intercept values along with "constrained slope" - 182 // the slope when the y-intercept is forced to zero. Optionally, apply the fit to the original data values and computed the 183 // standard deviation (Sigma) of the resulting drift-removed dataset. Drift-removed data values are discarded, original data elements are unmodified 184 // Example 1: do a linear fit during calibration to compute an angle - "Sigma" is not needed 185 // Example 2: do a linear fit on Dec values during a GA run - use the slope to compute a polar alignment error, use Sigma to estimate seeing of drift-corrected Dec values 186 // Returns a coefficient of determination, R-Squared, a form of correlation assessment 187 double GetLinearFitResults(double* Slope, double* Intercept, double* Sigma = NULL) const; 188 189 }; 190 191 class WindowedAxisStats : public AxisStats 192 { 193 bool autoWindowing = false; 194 int windowSize = 0; 195 void AdjustMinMaxValues(); 196 197 public: WindowedAxisStats()198 WindowedAxisStats() {}; 199 WindowedAxisStats(int AutoWindowSize); 200 ~WindowedAxisStats(); 201 202 // Change the window size of an active dataset - all stats will be adjusted accordingly to reflect the most recent <NewSize> elements 203 bool ChangeWindowSize(unsigned int NewWSize); 204 void RemoveOldestEntry(); 205 void AddGuideInfo(double DeltaT, double StarPos, double GuideAmt); 206 }; 207 208 #endif 209