I stumbled over this[1] page recently, which has a method that's apparently is better than the "traditional" Savitzky-Golay filters.
The idea seems to be to start with the desired frequency response, with lower frequencies close to the ideal differentiator, and higher frequencies tending smoothly to zero. This is then used to derive the filter coefficients through a set of equations.
The author generalizes it to irregularly sampled data near the end, so would be interesting to compare the approaches.
Just thought it'd throw it out there.
[1]: http://www.holoborodko.com/pavel/numerical-methods/numerical...
[1] https://www.sciencedirect.com/science/article/abs/pii/002192...
[2] https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=924...
I would've needed this recently for some data analysis, to estimate the mass of an object based on position measurments. I tried calculating the 2nd derivative with a Savitzky-Golay filter, but still had some problems and ended up using a different approach (also using a Kalman filter, but with a physics-based model of my setup).
My main problem was that I had repeated values in my measurements (sensor had a lower, non-integer divisible sampling rate than the acquisition pipeline). This especially made clear that np.gradient wasn't suitable, because it resulted in erratic switches between zero and the calculated derivative. Applying, np.gradient twice made the data look like random noise.
I will try using this library, when I next get the chance.
A common way to do that is to use a Savitzky-Golay filter [0], which does a similar thing - it can smooth out data and also provide smooth derivatives of the input data. It looks like this post's technique can also do that, so maybe it'd be useful for my ice strain-rate field project.
[0] - https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter
The usual implementation is as a convolution filter based on the premise that the underlying data is regularly sampled.
The pain in the arse occassional reality is missing data and|or present but glitched|spiked data .. both of which require a "sensible infill" to continue with a convolution.
This is a nice implementation and a potentially useful bit of kit- the elephant in the room (from my PoV) is "how come the application domain is irregularly sampled data"?
Generally (engineering, geophysics, etc) great lengths are taken to clock data samples like a metronome (in time and|or space (as required most)).
I'm assuming that your gridded GeoTIFF data field is regularly sampled in both the X and Y axis?
We had, in our geophysics application, a "pre computed" coefficient cache - the primary filters (central symmetric smoothing at various lengths) were common choices and almost always there to grab - missing values were either cheaply "faked" for Quick'NDirty displays or infilled by prediction filters that were S-G's computed to use existing points within the range to replace the missing value, that was either a look up from indexed filter cache or a fresh filter generation to use and stash in cache.
It's a complication (in the mechanical watch sense) to add, but with code to generate coefficients already existing it's really just looking at the generation times versus the hassle of indexing and storing them as created and the frequency of reuse of "uncommon" patterns.
Imagine you have a speed sensor eg. on your car and you would like to calculate the jerk (2nd derivative of speed) of your motion (useful in a range of driving comfort metrics etc.). The speed sensor on your car is probably not all that accurate, it will give some slightly randomly wrong output and it may not give that output at exactly 10 times per second, you will have some jitter in the rate you receive data. If you naiively attempt to calculate jerk by doing central differences on the signal twice (using np.gradient twice) you will amplify the noise in the signal and end up with something that looks totally wrong which you will then have to post process and maybe resample to get it at the rate that you want. If instead of np.gradient you use kalmangrad.grad you will get a nice smooth jerk signal (and a fixed up speed signal too). There are many ways to do this kind of thing, but I personally like this one as its fast, can be run online, and if you want you can get uncertainties in your derivatives too :)
The applications in this case involve self-driving cars, rocketry, homeostatic medical devices like insulin pumps, cruise control, HVAC controllers, life support systems, satellites, and other scenarios.
This is mainly due to a type of controller called the PID controller which involves a feedback loop and is self-balancing. The purpose of a PID controller is to induce a target value of a measurement in a system by adjusting the system’s inputs, at least some of which are outputs of the said controller. Particularly, the derivative term of a PID controller involves a first order derivative. The smoother its values are over time, the better such a controller performs. A problem where derivative values are not smooth or second degree derivative is not continuous, is called a “derivative kick”.
The people building these controllers have long sought after algorithms that produce at least a good approximation of a measurement from a noisy sensor. A good approximation of derivatives is the next level, a bit harder, and overall good approximations of the derivative are a relatively recent development.
There is a lot of business here. For example, Abbott Laboratories and Dexcom are building continuous blood glucose monitors that use a small transdermal sensor to sense someone’s blood glucose. This is tremendously important for management of people’s diabetes. And yet algorithms like what GP presents are some of the biggest hurdles. The sensors are small and ridiculously susceptible to noise. Yet it is safety-critical that the data they produce is reliable and up to date (can’t use historical smoothing) because devices like insulin pumps can consume it at real time. I won’t go into this in further detail, but incorrect data can and has killed patients. So a good algorithm for cleaning up this noisy sensor data is both a serious matter and challenging.
The same can be said about self-driving cars - real-time data from noisy sensors must be fed into various systems, some using PID controllers. These systems are often safety-critical and can kill people in a garbage in-garbage out scenario.
There are about a million applications to this algorithm. It is likely an improvement on at least some previous implementations in the aforementioned fields. Of course, these algorithms also often don’t handle certain edge cases well. It’s an ongoing area of research.
In short — take any important and technically advanced sensor-controller system. There’s a good chance it benefits from advancements like what GP posted.
P.S. it’s more solved with uniformly sampled data (i.e. every N seconds) than non-uniformly sampled data (i.e. as available). So once again, what GP posted is really useful.
I think they could get a job at pretty big medical and automotive industry companies with this, it is “the sauce”. If they weren’t already working for a research group of a self-driving car company, that is ;)
Give people clay/parchment/paper, however, and it becomes much easier to reliably evaluate an I term.
Example: https://xkcd.com/1205/ ; maybe each single time you do the task it seems like sanding out the glitches would be more trouble than it's worth, but a little record keeping allows one to see when a slight itch becomes frequent enough to be worth addressing. (conversely, it may be tempting to automate everything, but a little record keeping allows one to see if it'd obviously be rabbit holing)
One might say a part of the reason they have failed is because nature and people don't much care for the I term. These systems have feedback loops for sudden events, but increase the temperature slowly enough and the proverbial frog boils.
There are very many undercurrents in our world we do not see. So much that even when we think we understand and can predict their effects, we almost never take into account the entire system.
(the notion that we tend to overly ascribe stability and reproducibility to a system reminds me of Vilfredo Pareto having convinced himself that 80-20 was the invariable power law)
— We're taking turns to act as a sort of executive officer for the week
— Yes...
— But all the decisions of that officer have to be ratified at a special bi-weekly meeting
— Yes I see!
— By a simple majority, in the case of purely internal affairs
— Be quiet!
— But by a two-thirds majority, in the case of more major
— Be quiet! I order you to be quiet!
(I wonder how the Pythons themselves handled group decisions?)
* but see expander networks for a possible alternative network structure
I made a simulation that forecasted a greenhouse's temperature and humidity to help me understand the idea. I began by going over the basics of Gaussians and normal distribution once more. After that, I used NumPy and SciPy to develop the Kalman Filter in Python. To represent the system, I defined noise matrices (Q and R), a state transition matrix (F), and a control input matrix (B).
What’s your opinion on ensemble KF? We’d like to use that for parameter estimations. I saw unscented in your bayesfilter, but not ensemble, so I’m curious. Thanks!