You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: content/tutorials/empirical_mode_decomposition/empirical_mode_decomposition.qmd
+21-17Lines changed: 21 additions & 17 deletions
Original file line number
Diff line number
Diff line change
@@ -16,7 +16,7 @@ If signals oscillate regularly, Fourier methods are applicable.
16
16
These methods are well known and readily applied.
17
17
18
18
19
-
However, these methods fail to capture the complex characteristics of a signal.
19
+
However, these methods occasionally fail to capture the complex characteristics of a signal.
20
20
Irregular changes, varying amplitudes, and changing frequencies can result in derived measures which are no good representation of the data.
21
21
22
22
@@ -78,7 +78,7 @@ import matplotlib.pyplot as PLT
78
78
79
79
## gape angle data
80
80
81
-
For exploration, take this brief episode of beak gape angle of a canary cracking hemp seeds, courtesy of [Maja Mielke](https://orcid.org/0000-0001-6328-0589) ([*website*](http://mielke-bio.info/maja/blog)).
81
+
For exploration, take this brief episode of beak gape angle of a canary cracking hemp seeds, courtesy of [Maja Mielke](https://orcid.org/0000-0001-6328-0589) ([*website*](https://maja.mielke.ws/blog/)).
82
82
83
83
84
84
::: {.panel-tabset group="language"}
@@ -144,7 +144,7 @@ ShowPlot()
144
144
145
145
:::
146
146
147
-
This bird was in ["positioning" and "biting" phase](http://mielke-bio.info/maja/blog/01_eating_or_being_eaten), placing a hemp seed in the right position for cracking with the help of their upper beak, lower beak and tongue.
147
+
This bird was in ["positioning" and "biting" phase](https://doi.org/10.1242/jeb.244360), placing a hemp seed in the right position for cracking with the help of their upper beak, lower beak and tongue.
148
148
On the x axis of [@fig-data-py]/[@fig-data-r], you see the time (available in video frames, or normalized).
149
149
On the y axis, gape angle indicates (approximately) the angle between the beak tips and corner[^1].
150
150
@@ -160,7 +160,7 @@ You can clearly see an initial episode with large gape angle oscillations, and a
160
160
Engineers often talk about "peak amplitude" (as in: amplitude along the peaks) and quantify a **momentary amplitude** of a noticably oscillating signal ([see here](https://www.keysight.com/used/be/en/knowledge/guides/how-to-measure-amplitude-engineers-guide)).
161
161
162
162
This goes by the framework of ["Hilbert Transform"](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html#scipy.signal.hilbert) and ["Analytical Signal"](https://en.wikipedia.org/wiki/Analytic_signal).
163
-
All fascinating, yet we will here stick with the practical implementation.
163
+
All fascinating, yet we will here stick to the practical implementation.
164
164
165
165
166
166
To find the peak amplitude of the signal, one obviously has to detect the **peaks** (i.e. local maxima) of the signal.
@@ -295,7 +295,7 @@ One which is used by the reference implementation is `pchip`, "Piecewise Cubic"[
295
295
::: {.panel-tabset group="language"}
296
296
### R
297
297
298
-
In R, the `pracma` and `signal` libraries contain interpoltion functions which are not accurate (copies of matlab).
298
+
In R, the `pracma` and `signal` libraries contain interpolation functions which are not accurate (copies of matlab).
299
299
300
300
Then, there is the `interpolators` library, yet for some weird design choice it does not offer the choice to extrapolate beyond the data range (though splines can do that easily, see Python).
301
301
@@ -596,7 +596,7 @@ These functions can help a lot to play around with the peak finding parameters a
596
596
# Example II: Random Walk on Water
597
597
## anachronistic measures
598
598
599
-
Ultimately, I would like to apply these functions to quasi-continuous water level measurements from `watina`.
599
+
Ultimately, I would like to apply these functions to quasi-continuous water level measurements from [`watina`](https://www.vlaanderen.be/inbo/datasets/watina/).
600
600
My goal is to slightly improve water level analysis, for the following reason.
601
601
602
602
@@ -620,7 +620,7 @@ the mean of the three lowest/highest ground water levels measured in bi-weekly m
620
620
621
621
622
622
The anachronism is that we still calculate `xG3` like this, despite the availability of high frequency sampled water levels.
623
-
We artificially pick bi-weekly values.
623
+
We artificially pick bi-weekly values, at 8:00 AM or at noon.
624
624
We disregard the continuous, temporally periodic nature of the phenomenon.
625
625
We choose an arbitrary sampling cadence.
626
626
We ignore measurement uncertainty.
@@ -1059,7 +1059,7 @@ PLT.show();
1059
1059
As with the random walks above, unguided peak detection *should* find each tiny local peak, therefore initially extracting white noise where it is present.
1060
1060
1061
1061
1062
-
## Guided Mode Search
1062
+
## Guided Mode Search (*advanced*)
1063
1063
1064
1064
Experimenting with the peak detection parameters is worth a try.
1065
1065
(Because this works much better in Python, I will omit the R code here.)
@@ -1069,6 +1069,7 @@ First, search the long-term oscillations, by setting `prominence`, `width`, or `
1069
1069
In this special case, `width` will exclude the local minima: the lower peaks seem to be rather narrow on this water hole.
1070
1070
If this is an issue, remember that you can control peak finding for maxima and minima separately.
1071
1071
However, a combination of distance (more than half a year; remember that peaks and troughs are found separately) and prominence (*meaningful* peak) did the job to find the yearly baseline and straighten out the data (which might not be what we want).
1072
+
This **effectively eliminates drift**.
1072
1073
1073
1074
1074
1075
```{python py-emd-wata-topdown1}
@@ -1121,8 +1122,10 @@ There are some obvious problems in this:
1121
1122
- wet summer ~2015: pronounced "dry" minimum was lacking; water levels around that year are dragged towards the zero
1122
1123
- peak chopping: some peaks, e.g. mid-2010 and 2016, are lost by smoothing
1123
1124
1125
+
Note that, in this particular example, drainage might be the cause of the asymmetrical (capped) distribution of water level values; EMD still handles it relatively well.
1124
1126
1125
-
Nevertheless, there is another valuable outcome of EMD: we get the lower and upper **envelope**!
1127
+
1128
+
There is another valuable outcome of EMD: we get the lower and upper **envelope**!
1126
1129
1127
1130
```{python py-wateremd-envelope}
1128
1131
#| label: fig-water-envelope-py
@@ -1140,10 +1143,11 @@ PLT.show();
1140
1143
1141
1144
```
1142
1145
1146
+
This envelope, if averaged over time, is somewhat related to the GxG's.
1143
1147
1144
1148
1145
1149
1146
-
Note that there is no guarantee that the same guiding parameters will work on every observation in your data set.
1150
+
There is no guarantee that the same guiding parameters will work on every observation in your data set.
1147
1151
In my experience, a "guided" (top-down or bottom-up) EMD approach requires a lot of fiddling with the parameters, possibly even case distinction.
1148
1152
1149
1153
@@ -1162,13 +1166,13 @@ PLT.show();
1162
1166
1163
1167
```
1164
1168
1165
-
However, these peak detection controls are unavailable in R, so we might better stick with the defaults anyways (which usually work).
1169
+
The peak detection controls are unavailable in R (at least I did not find them), so we might better stick with the defaults anyways (which usually work).
1166
1170
1167
1171
1168
1172
## Summary: Water Level EMD
1169
1173
1170
1174
::: {.callout-note}
1171
-
To "wrap up" (envelope-pun), some observations:
1175
+
To "wrap up" (envelope-pun intended), some observations:
1172
1176
1173
1177
- "guided EMD": yearly oscillations are found first by selecting for peak characteristics
1174
1178
- noise can be extracted either before or after; EMD is one way of naturally smoothing a signal
@@ -1181,17 +1185,17 @@ And the most important take-home message:
1181
1185
:::
1182
1186
1183
1187
1184
-
One goal of my application of EMD was to find a better way to extract a yearly range of water levels, to replace the anachronistic `LG3` and `HG3` calculations.
1188
+
One goal of my application of EMD was to find a better way to extract a yearly range of water levels, to replace the anachronistic `LG3` and `HG3` calculations (which, nevertheless, worked surprisingly well on random walks).
1185
1189
With plain application of EMD (i.e. no pre-processing of the data), the envelope seems to be a promising aspect for further inspection.
1186
1190
1187
-
The reason that water levels turned out to be non-ideal for EMD application is that they lack regular oscillations.
1188
-
They are not symmetric (winter wet plateau; summer dry dip), and not necessarily regular (e.g. "wet summer").
1189
-
Generally, water level measurements such as the first example above might be more usefully approached with [wavelets](https://docs.scipy.org/doc/scipy-1.12.0/reference/signal.html#wavelets) to find the summer minima.
1191
+
The cases in which water levels turned out to be non-ideal for EMD application is when they lack regular oscillations.
1192
+
They are not symmetric (e.g. winter wet plateau due to drainage, summer dry dip, drift, ...), and not necessarily regular (e.g. "wet summer").
1193
+
Water level measurements such as the first example above might be more usefully approached with [wavelets](https://docs.scipy.org/doc/scipy-1.12.0/reference/signal.html#wavelets) to find the summer minima.
1190
1194
This is just my [almost-uneducated guess](http://mielke-bio.info/falk/posts/27.cycle_extraction/#orgac5a9c6), based on the visual form of the curves.
1191
1195
CWT (Continuous Wavelet Transform, [e.g. in R](https://www.rdocumentation.org/packages/Rwave/versions/2.6-5/topics/cwt)) is a great subject for another tutorial.
1192
1196
1193
1197
1194
-
As a reminder, there are many tools to consider for signal analysis.
1198
+
As a general reminder, there are many tools to consider for signal analysis.
1195
1199
I hope this tutorial could help to bring EMD to your personal repertoire.
0 commit comments