Interactive companion for fitting splines to load trajectories, projecting samples, and reading windowed Fourier transforms and wavelet scalograms
A load carried up a curved tree leaves a trajectory whose shape is dominated by the tree itself. To detect anything behavioral riding on top of that shape – pre-fork hesitation, side-to-side deliberation, drift – the tree’s own contribution has to be separated out first. This widget works through three ways of doing that, each giving up something the one before it relied on.
d(s), with the smoothing chosen either from what you already know about the
tree or from patterns in the data.x(z) against height calls the low spatial frequencies tree and the rest ants,
and never estimates where the tree is.The trajectory is modeled as an autocorrelated random walk along the centerline, clipped to the trunk surface – every sample sits on the tree, mimicking samples from a single video of a load being cooperatively transported. Approaching the fork, the walk becomes both wider and quicker; the change in speed, not the change in width, is what a spatial-frequency cutoff can separate from the tree.
t is the normalized trajectory parameter:
t = 0 at the base of the trunk, t = 1 at the
tip of the right branch, and t = 0.5 at the fork.
The lateral coordinate follows an AR(1) random walk, with an elevated wobble in a window
around t = 0.5 (the “deliberation” signal).
Green curve is the cubic B-spline fitted to the samples via least squares.
Red lines are perpendicular deviations di from each sample to its closest point on the spline.
What is DF? These points are fit with a cubic spline, a smooth chain of cubic polynomials, and DF counts the degrees of freedom, meaning the free parameters available to shape that fit. The Choosing spline parameters tab takes the idea apart in full, and the Perpendicular projection tab shows what the residuals off such a fit are measuring. In short, a single cubic (a polynomial reaching an exponent of 3) carries 4 free parameters, the coefficients of 1, t, t², and t³, and so DF = 4 leaves zero interior knots, one unbroken polynomial spanning the whole trajectory. The number of interior knots is DF − 4, that is DF − (degree + 1): DF = 5 has 1 knot, DF = 6 has 2, and DF = 14 has 10. Toggle Spline knots to see them marked.
di against arc length si
of the foot on the spline. For a Y-tree, the vertical orange line marks the arc-length location of the
tree's fork projected onto the fitted spline. The tree's fork is at a fixed geometric
point, but its arc-length coordinate along the spline shifts slightly with each refit because
the spline itself shifts. Dashed gray lines indicate ±1σ
of the residual estimated away from the deliberation window.
If the spline captures the tree's shape – the smooth path the trajectory would
follow without behavior – the residuals d(s) contain everything else:
autocorrelated noise, drift, hesitation. Pre-fork wobble shows up as elevated residual variance
around the fork.
The smoothing is a Goldilocks choice: too few degrees of freedom and the spline cuts across the bend, dumping tree shape into the residual; too many and it threads through the deliberation, hiding the signal. See the Data-driven parameters tab.
On the smoothing parameter. How closely the spline is allowed to follow the trajectory
rests on a single dial. Here you set it directly, as the DF of the fit, which the
Choosing spline parameters tab takes apart. Some packages
(mgcv::gam, say) instead have you set how heavily wiggliness is charged for and
let the flexibility follow from that; the Data-driven parameters tab shows why those are
the same choice approached from opposite ends.
library(splines) # Inputs # ants: one row per video frame, columns t (trajectory parameter running # 0 to 1), x (lateral position), z (height up the tree). # Fit cubic B-splines to x and z independently as functions of t. df <- 6 fit_x <- lm(x ~ bs(t, df = df, degree = 3), data = ants) fit_z <- lm(z ~ bs(t, df = df, degree = 3), data = ants) # Dense evaluation of the curve r(t) = (x_hat, y_hat) t_grid <- seq(min(ants$t), max(ants$t), length.out = 600) curve <- data.frame( t = t_grid, x = predict(fit_x, newdata = data.frame(t = t_grid)), z = predict(fit_z, newdata = data.frame(t = t_grid)) ) # Cumulative arc length along the curve ds <- sqrt(diff(curve$x)^2 + diff(curve$z)^2) curve$s <- c(0, cumsum(ds))
library(splines)
library(dplyr)
library(tibble)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
df <- 6
fit_x <- lm(x ~ bs(t, df = df, degree = 3), data = ants)
fit_z <- lm(z ~ bs(t, df = df, degree = 3), data = ants)
curve <- tibble(t = seq(min(ants$t), max(ants$t), length.out = 600)) |>
mutate(
x = predict(fit_x, newdata = pick(t)),
z = predict(fit_z, newdata = pick(t))
) |>
mutate(
dx = x - lag(x, default = first(x)),
dz = z - lag(z, default = first(z)),
s = cumsum(sqrt(dx^2 + dz^2))
) |>
select(-dx, -dz)
Each sample pi = (xi, zi) is mapped to two scalars:
a signed deviation di from the fitted curve and an arc length
si giving its position along the curve.
1. Find the foot. The nearest point on the curve is
ti∗ = argmint∈[0,1] ‖ pi − r(t) ‖where ‖ · ‖ is the Euclidean distance – the ordinary straight-line distance √[(x1−x2)² + (z1−z2)²].
In practice we evaluate r(t) at a few hundred equally spaced parameter
values and pick the index where the squared distance is minimized.
2. Compute the signed scalar distance. Let τ(ti∗) be the unit tangent at the foot and θi the angle from τ(ti∗) to (pi − r(ti∗)), measured counterclockwise as shown in the figure. Define the sign factor sgn(θi) by embedding both vectors in 3D (zero z-components) and extracting the z-component of their cross product:
sgn(θi) = (0,0,1) · τ(ti∗) × (pi − r(ti∗)) ‖τ(ti∗)‖ · ‖pi − r(ti∗)‖Because r(ti∗) is the nearest point on the curve to pi, the displacement (pi − r(ti∗)) is perpendicular to τ(ti∗), and so sgn(θi) = ±1 exactly. The signed scalar distance is then:
di = ‖pi − r(ti∗)‖ · sgn(θi)Positive on the left of the oriented curve, negative on the right – though this sign convention is arbitrary. Reversing the cross-product order to (pi − r(ti∗)) × τ(ti∗) flips the sign and labels the right side positive instead.
3. Tag with arc length. The companion coordinate is
si = ∫0ti∗ ‖ r′(u) ‖ du
computed numerically as a cumulative sum of chord lengths over the same dense sequence of parameter values.
Plotting di against si gives the residual signal
used in the first tab.
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
# curve: the dense spline evaluation built by the Detrending overview
# snippet, columns x, z, s (cumulative arc length). This snippet
# continues that one rather than standing alone.
# Project samples (x, z) onto a dense evaluation of the fitted
# spline (curve$x, curve$z) and compute the SCALAR d_i and s_i.
K <- nrow(curve)
# Unit tangent tau at each curve point (central difference)
i_prev <- pmax(seq_len(K) - 1, 1)
i_next <- pmin(seq_len(K) + 1, K)
tx <- curve$x[i_next] - curve$x[i_prev]
tz <- curve$z[i_next] - curve$z[i_prev]
mag <- sqrt(tx^2 + tz^2)
tx <- tx / mag; tz <- tz / mag
# For each sample, find nearest curve point (foot)
project_one <- function(px, pz) which.min((curve$x - px)^2 + (curve$z - pz)^2)
i_foot <- mapply(project_one, ants$x, ants$z)
# Signed scalar distance: 2D cross product of unit tangent with gap vector.
# Equals ||p_i - r(t_i*)|| * sin(theta_i) where theta_i is the signed angle
# from tau to the displacement. Positive on the left of the oriented curve.
d_i <- tx[i_foot] * (ants$z - curve$z[i_foot]) -
tz[i_foot] * (ants$x - curve$x[i_foot])
s_i <- curve$s[i_foot]
resid <- data.frame(s = s_i, d = d_i)library(dplyr)
library(purrr)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
# curve: the dense spline evaluation built by the Detrending overview
# snippet, columns x, z, s (cumulative arc length). This snippet
# continues that one rather than standing alone.
# Add unit tangent columns to the curve tibble
curve <- curve |>
mutate(
tx = lead(x, default = last(x)) - lag(x, default = first(x)),
tz = lead(z, default = last(z)) - lag(z, default = first(z)),
mag = sqrt(tx^2 + tz^2),
tx = tx / mag,
tz = tz / mag
)
# Project each sample; d_i = ||p_i - r(t_i*)|| * sin(theta_i)
# computed directly as the 2D cross product tau x delta.
resid <- ants |>
mutate(
i_foot = map2_int(x, z, \(px, pz) which.min((curve$x - px)^2 + (curve$z - pz)^2)),
d = curve$tx[i_foot] * (z - curve$z[i_foot]) -
curve$tz[i_foot] * (x - curve$x[i_foot]),
s = curve$s[i_foot]
)When the relationship between two quantities does not appear to follow a straight line, the usual first move is to model it with a polynomial: a weighted sum of powers of whatever you are treating as the input. That is a modeling choice and not a claim about what the relationship really is; a polynomial is simply a shape flexible enough to be bent into agreement with the data. One power alone gives a straight line, adding a squared term lets the curve bend, adding a cubed term lets that bend reverse, and the weights are what pick one shape out of the family.
A cubic is the polynomial that stops at the third power, and so four weights describe it
completely. Throughout this widget z is the height up the tree and x
is the load’s lateral position, and it is the lateral position we are modeling as the
height changes, and so a cubic here reads:
x = a + b z + c z2 + d z3
There is nothing else to adjust. a shifts the whole curve sideways, b
tilts it, and c and d bend it; the four panels below are that one
expression under four different sets of numbers. What no choice of the four can do is wiggle
much: a cubic turns at most twice, and reverses which way it curves at most once.
Four cubics, each labeled with the numbers that produce it. Horizontal is the height
z; vertical is the lateral position x; the faint line is
x = 0. The second has no z³ term at all, which is
what leaves it a plain parabola – the cubic family contains the simpler shapes too.
A load carried up a tree wanders far more than that, and so one cubic will not do. Cut the height range into stretches and give each stretch its own cubic. The cut points are the knots, and they are the real decision: they say where the curve is allowed to change what it is doing.
Fitted independently, neighboring cubics would disagree wherever they meet – a step, a kink, or an abrupt change of curvature. A spline is that same chain with the disagreements forbidden: at every join the two cubics must match in height, in slope, and in curvature. Three conditions per join, and that is the whole of the machinery.
Which is where the parameter count comes from. With n knots there are
n+1 stretches, each carrying its own four numbers, so 4(n+1) to begin
with; the three conditions at each of the n joins take 3n away:
DF = 4(n+1) − 3n = n + 4
So the fit has n + 4 numbers free to move, its degrees of freedom. The 4 is
the four numbers of a single cubic, sitting there because you always have at least one; each
knot after that buys exactly one more.
Why the knots are a filtering choice. Each stretch is one cubic, and so it can turn once. Reproducing a wiggle therefore takes about two stretches, one for each way it turns. Wiggles longer than that come through the fit; shorter ones get flattened, and the second panel below lets you watch the changeover. That is the whole heuristic: choosing knots per meter is choosing how much detail per meter you are willing to treat as trend rather than signal.
library(splines)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
# df = interior knots + 4 for a cubic, so df = 6 puts 2 knots inside the range.
fit <- lm(x ~ bs(z, df = 6, degree = 3, intercept = TRUE), data = ants)
ants$trend <- fitted(fit)
ants$resid <- resid(fit)
# Left alone, bs() puts knots at quantiles of z, so their spacing follows the
# data. Pass them explicitly for uniform spacing, so the fit smooths away
# detail finer than one knot gap evenly over the whole tree.
n_int <- 2
brk <- seq(min(ants$z), max(ants$z), length.out = n_int + 2)
fit2 <- lm(x ~ bs(z, knots = brk[-c(1, length(brk))],
degree = 3, intercept = TRUE), data = ants)
plot(ants$z, ants$x, pch = 16, col = "gray70",
xlab = "height z", ylab = "lateral x")
lines(ants$z, fitted(fit2), col = "#6b21a8", lwd = 2)
abline(v = brk[-c(1, length(brk))], lty = 3, col = "gray50")library(dplyr)
library(broom)
library(ggplot2)
library(splines)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
n_int <- 2
brk <- seq(min(ants$z), max(ants$z), length.out = n_int + 2)
knots <- brk[-c(1, length(brk))]
fitted_ants <- lm(x ~ bs(z, knots = knots, degree = 3, intercept = TRUE),
data = ants) |>
augment(data = ants) |>
rename(trend = .fitted, resid = .resid)
ggplot(fitted_ants, aes(z, x)) +
geom_point(color = "gray70") +
geom_line(aes(y = trend), color = "#6b21a8", linewidth = 1) +
geom_vline(xintercept = knots, linetype = 3, color = "gray50") +
labs(x = "height z", y = "lateral x") +
theme_minimal()Where DF can come from. The previous tab set it from what you already know about the tree: a bend scale, or a knot density per meter coarse enough that anything finer must be the ants. That is the right route whenever such prior knowledge exists. When it does not, the alternative is to let the trajectory itself decide, reading DF off patterns in the data. That is model selection, and it comes in two families, differing in whether DF is set directly or allowed to emerge.
Let the penalty settle it. A penalized spline (e.g. mgcv::gam) keeps
plenty of knots throughout and charges the fit for how much it wiggles, letting the size of
that charge decide what it can afford: a heavy charge buys a smooth curve with few effective
degrees of freedom, a light one buys a flexible curve with many. The charge is usually written
λ, and raising DF is the same move as easing λ off, and so the two are
parameterizations of one trade-off rather than rival methods. GCV then picks λ
without being told, by minimizing an approximation of leave-one-out prediction error, and the
effective DF settles only as high as the data warrant. Nothing has to be swept, which is why
this is the shorter road when you have no prior knowledge to bring.
Or sweep DF and score every fit. Ordinary least squares, used here and in
splines::bs() + lm(), fixes DF before it fits, and so
choosing DF means fitting at every DF on a grid and scoring the results. That needs a score,
and the rest of this tab builds one. It is more setup than the penalized route, and it earns
that setup by making the failure visible: you can watch overfitting begin.
Define the residual sum of squares between two parametric curves sampled at N common parameter values as
RSS(A, B) ≜ ∑k=1N ‖ Ak − Bk ‖2Both quantities in the plot below compare the fitted spline to a reference set of N points evaluated at the same parameter values:
• RSS(spline, data) — spline points vs. sample positions. Decreases monotonically with DF (more flexibility always reduces in-sample error, reaching zero when DF ≥ N). It has no minimum, and so it cannot choose DF on its own. RSS(spline, data) is shown not because it picks DF, but as a foil: its relentless decrease illustrates why a monotone fit criterion is insufficient and a complexity penalty is needed.
• RSS(spline, tree) — spline points vs. true tree-centerline positions (oracle). For a Y-shaped tree it is U-shaped: at low DF the spline cannot follow the fork geometry (underfitting); at high DF it chases trajectory noise away from the true shape (overfitting), and the minimum df* identifies the ideal smoothing level. For a straight (I-shaped) tree the geometry is linear – fully captured at df = 4 – so there is no underfitting regime and RSS(spline, tree) rises monotonically: higher DF only lets the spline drift further from the straight backbone.
How to use both curves together. The two curves diverge as DF grows: RSS(spline, data) keeps falling whereas RSS(spline, tree) passes its minimum and starts to rise. The onset of this divergence marks the overfitting threshold – at that DF, the spline begins fitting behavioral noise rather than tree geometry. Practically, choose DF in the neighborhood of df* (or slightly to its left for a conservative smoothing). The widening gap between the two curves tells you how much of the fit improvement is noise-fitting rather than tree-shape recovery: a narrow gap means most of the fit gain is genuinely structural; a wide gap means the spline is mostly chasing noise.
Can you get the tree skeleton from real video? Yes – the trunk centerline can be extracted by image segmentation or frame-to-tree registration. When that is possible, RSS(spline, tree) is directly usable, not just an oracle. When it is not available, AIC, BIC, and GCV replace the oracle: instead of measuring how far the spline is from the true skeleton, they penalize the in-sample RSS by the number of parameters, discouraging the spline from chasing noise.
• AIC = N log(RSS/N) + 2 DF
(Akaike – lighter penalty, slightly favors richer models).
• BIC = N log(RSS/N) + log(N) DF
(heavier penalty, sparser; usually preferable for large N).
Both AIC and BIC are minimized over a grid of DF values.
AIC, BIC, and GCV all penalize complexity for the same reason: a spline with enough DF can pass through every sample (RSS = 0), but that spline is useless – it captures noise, not tree shape. The penalty mimics the U-shape of RSS(spline, tree) that arises for non-straight trees, without needing the true skeleton.
df overlaid on the same trajectory. The
current df (from Tab 1's slider) is drawn in blue; the fixed reference values
are gray. The dark-brown centerline running through the tree is the true skeleton backbone
used as ground truth for RSS(spline, tree).
df; the dotted green line
marks df* = argmin RSS(spline, tree). Hit
Average over 30 seeds to smooth the noise realizations and see the underlying
U-shape clearly.
library(splines) # provides bs()
library(mgcv)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
# skel: the tree centerline at each sample's t, columns x, z. Simulation
# only, as real data has no known centerline; Option 1 below is the
# only part that needs it.
# bs(t, df = k, degree = 3) builds a cubic B-spline design matrix: an n x k
# matrix whose columns are B-spline basis functions evaluated at each t value.
# Passing it to lm() fits by OLS over that k-dimensional function space —
# the same least-squares cubic B-spline with k free parameters used throughout
# this widget. Interior knots are placed at equally spaced quantiles of t;
# df = k = (degree + 1) + number_of_interior_knots = 4 + interior knots.
# Option 1: compare RSS(spline, data) vs RSS(spline, tree) across a DF grid.
# `skel` holds the tree-centerline (x, z) at each sample's t value.
# RSS(spline, data) always decreases; RSS(spline, tree) is U-shaped for Y-trees.
rss_df <- function(df) {
fit_x <- lm(x ~ bs(t, df = df, degree = 3), data = ants)
fit_z <- lm(z ~ bs(t, df = df, degree = 3), data = ants)
xh <- predict(fit_x, newdata = ants)
zh <- predict(fit_z, newdata = ants)
c(data = sum((ants$x - xh)^2 + (ants$z - zh)^2),
tree = sum((skel$x - xh)^2 + (skel$z - zh)^2))
}
rss_mat <- sapply(4:20, rss_df)
colnames(rss_mat) <- 4:20
df_star <- as.integer(colnames(rss_mat)[which.min(rss_mat["tree", ])])
# Option 2a: sweep DF, pick by AIC or BIC
ic <- sapply(4:20, function(df) {
fx <- lm(x ~ bs(t, df = df, degree = 3), data = ants)
fz <- lm(z ~ bs(t, df = df, degree = 3), data = ants)
c(aic = AIC(fx) + AIC(fz),
bic = BIC(fx) + BIC(fz))
})
colnames(ic) <- 4:20
df_aic <- as.integer(colnames(ic)[which.min(ic["aic", ])])
df_bic <- as.integer(colnames(ic)[which.min(ic["bic", ])])
# Option 2b: penalized cubic regression spline, smoothness by GCV
gx <- gam(x ~ s(t, bs = "cr"), data = ants, method = "GCV.Cp")
gz <- gam(z ~ s(t, bs = "cr"), data = ants, method = "GCV.Cp")
edf <- gx$edf + gz$edf # effective DF from GCVlibrary(splines) # provides bs()
library(mgcv)
library(dplyr)
library(purrr)
library(tibble)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
# skel: the tree centerline at each sample's t, columns x, z. Simulation
# only, as real data has no known centerline; Option 1 below is the
# only part that needs it.
# bs(t, df = k, degree = 3) constructs a cubic B-spline design matrix: k basis
# functions evaluated at each t, used by lm() as predictors for OLS fitting.
# This produces the same least-squares cubic B-spline (k free parameters) as
# the widget's interactive spline. Interior knots fall at t quantiles;
# df = k = 4 + number_of_interior_knots.
# Option 1: compare RSS(spline, data) vs RSS(spline, tree) across a DF grid.
rss_sweep <- tibble(df = 4:20) |>
mutate(
fit_x = map(df, \(k) lm(x ~ bs(t, df = k, degree = 3), data = ants)),
fit_z = map(df, \(k) lm(z ~ bs(t, df = k, degree = 3), data = ants)),
rss_data = map2_dbl(fit_x, fit_z, \(fx, fz) {
xh <- predict(fx, newdata = ants)
zh <- predict(fz, newdata = ants)
sum((ants$x - xh)^2 + (ants$z - zh)^2)
}),
rss_tree = map2_dbl(fit_x, fit_z, \(fx, fz) {
xh <- predict(fx, newdata = ants)
zh <- predict(fz, newdata = ants)
sum((skel$x - xh)^2 + (skel$z - zh)^2)
})
)
df_star <- rss_sweep |> slice_min(rss_tree) |> pull(df)
# Option 2a: sweep DF, pick by AIC / BIC
ic <- rss_sweep |>
mutate(
aic = map2_dbl(fit_x, fit_z, \(fx, fz) AIC(fx) + AIC(fz)),
bic = map2_dbl(fit_x, fit_z, \(fx, fz) BIC(fx) + BIC(fz))
)
df_aic <- ic |> slice_min(aic) |> pull(df)
df_bic <- ic |> slice_min(bic) |> pull(df)
# Option 2b: GCV-chosen smoothness via mgcv::gam
gx <- gam(x ~ s(t, bs = "cr"), data = ants, method = "GCV.Cp")
gz <- gam(z ~ s(t, bs = "cr"), data = ants, method = "GCV.Cp")
edf <- gx$edf + gz$edf
Instead of fitting a spline and subtracting it, we can move directly into the spatial-frequency
domain. Let x(z) be the observed lateral coordinate as a function of height
y, resampled onto a uniform grid. Nothing about the tree is used: the tree’s
shape is precisely what the spline tabs have to estimate, and the point of this route is to
avoid estimating it. The trunk’s gentle curvature, the branch’s sideways ramp, and
the trajectory’s slow drift all contribute power at low spatial frequencies;
deliberation contributes faster side-to-side movement, and therefore power at higher
frequencies. A windowed Fourier transform of x(z) exposes both at once.
Nothing is subtracted off first. That is the whole point: the transform is what does the separating, and so the trend has to be left in for the divide to have anything to separate, and the tree shows up as the bright band along the left edge. The one exception is each window’s mean, which is not a trend at all, being wherever the horizontal origin happens to sit. Leaving it in would put a single DC bin about 17 times above anything past the divide and flatten the rest of the plot to black.
The slow trajectory shape has some characteristic length scale Lbend.
The divide at spatial frequency fcut = 1 / Lbend reaches the
same end as a spline detrend without fitting or subtracting anything: power below it is the
tree, power above it is candidate behavioral structure. Drag the Lbend
slider to move the divide. What you are looking for is energy crossing to the right of
the dashed line as the load approaches the fork.
Why slice by height, not distance traveled. Height is directly observable and, unlike the distance traveled along the trajectory, is unaffected by the wobble being measured. Parameterizing by path length would let a bout of wide side-to-side movement advance the axis faster, stretching that bout to a longer apparent wavelength and pushing it the wrong way across the cutoff. Height also lets this panel share a vertical axis with the tree beside it, and so the dashed fork line falls at the same place in both and any height can be read straight across.
Two things to be wary of. The fork is a kink in the trend itself, and so the tree alone puts
a broadband smudge at that height; its energy falls off like 1/f², and so most of
it should stay left of the cutoff, but the streak at the fork row is not evidence of behavior.
And the frequency axis stops at the Nyquist limit of the raw samples, not beyond it: resampling
onto a finer grid would only interpolate, and interpolation corners are broadband, producing a
floor that is invisible where the signal is strong and dominant wherever it is weak. That is how
a quiet, straight stretch of trunk can be painted as high frequency, especially under per-slice
normalization, which rescales an almost empty window to full brightness.
x(z) used
by the spectrogram on the right is the horizontal coordinate of these samples.
The dashed orange fork line is drawn at the same height in both panels, and on the
same pixel row wherever the two sit side by side.
x(z), as a function of height z. No knowledge of the tree is
used, and nothing is detrended: the trend stays in, which is why the bright band on the far
left is the tree itself. Only each window’s mean is subtracted, and that is an arbitrary
coordinate origin rather than a trend.
Vertical axis is height z, shared with the tree on the left; horizontal axis
is spatial frequency, stopping at the raw samples’ Nyquist limit. The horizontal orange
line marks the fork. The dashed white vertical line is the predicted divide
fcut = 1/Lbend in spatial frequency between power
contributions from the smoothly changing tree (left of the line) and from tighter changes
intrinsic to the ant dynamics (right of the line). The faint dashed lines marked
edge sit W/2 (half of the window length) from the top and bottom of the
tree, and more specifically of the trajectory data. Between the edges of the tree and those
edge lines, part of the W window is not reflective of the trajectory
data, and so read those bands with caution.
hop:
drop it to 1 and the trace fills in completely. That added density is cosmetic, as
neighboring windows then share all but one of their samples. What actually blurs the bulge
is W, every point being an average over a W-tall window, and no
hop setting recovers what a long window has already smeared. Raise
W and watch the bulge spread even at hop = 1, then
compare against the same panel on the next tab, which has neither knob to set.
Window length W — number of samples in each STFT slice. Larger W
gives better frequency resolution but coarser localization in z; a Heisenberg-style
trade-off, and one you can see directly, as the edge bands at the top and bottom of the
spectrogram are each W/2 deep and grow as you widen the window. The uniform grid
K is not a free parameter: it is set so the grid spacing matches the coarsest
spacing of the raw samples in z, which is what fixes the right-hand end of the
frequency axis. Both K and the resulting Nyquist limit are reported under the
Lbend slider, along with a note if W had to be capped
to fit the grid.
Hop — shift between consecutive windows. Smaller hop means more overlap (smoother
spectrogram, more compute). Default hop = 4 with
W = 20 gives 80% overlap. The default W is set so a window
spans about the same height as a deliberation bout: widen it much beyond that and the bout is
smeared across many slices, narrow it much below and too few frequency bins fall below
fcut to read the split.
Window type — Hann and Hamming taper smoothly to (near) zero at the edges, suppressing spectral leakage from window boundaries; their side-lobe behavior differs slightly. Rectangular applies no taper – sharp edges leak energy across frequencies. Gaussian tapers smoothly and is optimal in the joint time-frequency uncertainty sense (Gabor's bound).
Lbend — characteristic length scale of slow trajectory shape (the trunk's
bend, an ant's drift toward one side, anything you want to call “trend”).
The cutoff fcut = 1/Lbend partitions the spectrum:
below is trend, above is candidate signal.
log magnitude — switch from linear to log10(1+9v) scaling
of the spectrogram so weaker components become visible.
per-slice normalization — rescale each z-slice to its own max,
useful when overall power varies along the trunk. Read it with care: a slice whose absolute
power is near zero gets its noise floor promoted to full brightness, and so bright high-frequency
bands on quiet, straight stretches are the expected artifact of this option rather than a
finding. Cross-check anything interesting with the option off.
library(signal)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
# Lateral coordinate x resampled against HEIGHT z. Grid spacing matches the
# coarsest raw spacing in z, so the frequency axis stops where the samples stop
# supporting it; a finer grid would only interpolate, and interpolation corners
# are broadband.
K <- round(diff(range(ants$z)) / max(diff(ants$z))) + 1
z_grid <- seq(min(ants$z), max(ants$z), length.out = K)
xy <- approx(ants$z, ants$x, xout = z_grid)$y
dz <- diff(z_grid)[1]
W <- 20; hop <- 4 # W spans about one deliberation bout
sp <- specgram(xy, n = W, Fs = 1/dz, window = hanning(W), overlap = W - hop)
# specgram reports frame STARTS in sp$t, as offsets from the first sample.
# Shift to frame centers, expressed as height:
z_c <- min(ants$z) + sp$t + (W %/% 2 - 1) * dz
# Drop the DC row. It is the window mean, i.e. wherever the x-origin happens to
# sit, and it runs ~32x the rest of the spectrum. Dropping it does the same job
# as subtracting each window's mean: the next bin up is inflated only ~1.4x by
# leakage and the bins beyond that not at all. The trend itself is NOT removed,
# and shows up as the bright band at the left.
mag <- Mod(sp$S)[-1, ]
f <- sp$f[-1]
z_fork <- 1 # height of the geometric fork
L_bend <- 0.15 # length scale treated as tree rather than ant
f_cut <- 1 / L_bend # the divide, in 1/length
# VERTICAL = z, so it reads against a drawing of the tree
image(x = f, y = z_c, z = mag,
col = hcl.colors(64, "viridis"),
xlab = "spatial frequency (1/length)",
ylab = "height z")
abline(h = z_fork, col = "darkorange", lwd = 2)
abline(v = f_cut, col = "white", lty = 2, lwd = 1.5)
# This covers only the middle ~67% of the tree, and that is not a bug: specgram
# places a frame only where a whole window fits, so nothing is centered within
# W/2 of either end. Widen W and the covered band shrinks further. The widget
# instead holds the last observed value past each end, so the axis can be shared
# with a drawing of the tree, and marks those bands, as part of the window there
# is not real data. Zero-filling reads as the natural choice and is worse: the
# branch tip sits near x = 0.4, so zeros put a step there that lights up every
# frequency. Note that WaveletComp on the next tab DOES zero-pad internally,
# which is why its snippet pads by hand before calling it.library(signal)
library(tibble)
library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
# Uniform grid in HEIGHT, spaced to match the coarsest raw spacing so the
# frequency axis never runs past what the samples support.
K <- round(diff(range(ants$z)) / max(diff(ants$z))) + 1
z_grid <- seq(min(ants$z), max(ants$z), length.out = K)
xy <- approx(ants$z, ants$x, xout = z_grid)$y
dz <- diff(z_grid)[1]
W <- 20; hop <- 4 # W spans about one deliberation bout
win <- hanning(W)
pad <- W %/% 2
xp <- c(numeric(pad), xy, numeric(pad))
# Frames centered on grid points, each window's mean removed
spec_df <- tibble(ctr = seq(1, K, by = hop)) |>
mutate(
z = z_grid[ctr],
mag = map(ctr, \(i) {
seg <- xp[i:(i + W - 1)]
tibble(f = (0:(W %/% 2)) / (W * dz),
magnitude = Mod(fft((seg - mean(seg)) * win))[1:(W %/% 2 + 1)])
})
) |>
unnest(mag)
# Vertical axis = height z, horizontal axis = spatial frequency
ggplot(spec_df, aes(x = f, y = z, fill = magnitude)) +
geom_raster() +
geom_hline(yintercept = 1, color = "darkorange", linewidth = 1.2) +
geom_vline(xintercept = 1 / 0.15, color = "white", linetype = "dashed") +
scale_fill_viridis_c() +
labs(x = "spatial freq (1/length)", y = "height z") +
theme_minimal()
The previous tab has to commit to a window length W before it can look at
anything, and that one choice has to satisfy two incompatible demands. A window long enough to
contain a slow bend of the tree is longer than the deliberation bout you are trying to
localize, and so the bout gets smeared across many slices; a window short enough to pin the
bout down cannot represent the bend at all. There is no W that does both,
and no amount of extra data fixes it: both requirements are stated in height, and
shrinking the sample spacing changes neither.
A wavelet gets out of this by changing what the trajectory is compared against. Fourier analysis compares it against sine waves, and a sine wave runs forever in both directions. A tree of finite height cannot be compared against something infinite without first cutting a piece out of the data to hand it, which is all the window on the previous tab ever was, and the length of that piece is the choice that caused the trouble. A wavelet is instead a short packet: a few oscillations that rise, peak, and die away inside a finite stretch. It is already local, and so there is nothing to cut.
The Morlet wavelet used here is one popular choice among many, and a natural one to start with because it stays closest to Fourier: a sine wave multiplied by a Gaussian bump that fades it out at both ends. Widening or narrowing it stretches the oscillations and the bump together, in a single motion, and that is the difference that matters. On the previous tab, window length and frequency were two separate knobs, which is what forced the choice between them. Here they are tied to one number, the scale, and so a long packet is automatically a slow one and a short packet automatically a fast one. Coarse scales get the reach to see a slow bend, and fine scales stay sharply localized in height.
At every height, the stretched packet is lined up against the trajectory and compared with it, and what gets kept is how strongly the two resemble each other at that spot. Doing that at every height, for every scale in the family, gives one number per height per scale.
A scalogram is the picture that comes out, and it is the wavelet counterpart of the spectrogram on the previous tab: magnitude plotted against position (height, up the page) and scale (across the page, here converted to the equivalent spatial frequency so the two tabs read the same way). Where a spectrogram has one resolution everywhere, a scalogram trades resolution smoothly along its horizontal axis, fine in height at the right, fine in frequency at the left.
What is left to choose. Very little, which is the point. Position and scale are
both swept rather than set, as the transform visits every height at every rung of the
ladder. The width of the Gaussian envelope is not a separate setting either, being the
scale itself, and that is what ties the packet’s length to its frequency. One number
does remain, the Morlet w0 above, which fixes how many oscillations
fit under the envelope: raise it and each packet carries more cycles, which sharpens the
range of frequencies it responds to at the cost of blurring where in height it responds.
Lower it for the opposite trade. The value 6 is the usual default, and it is the one that
makes a packet’s scale and its wavelength very nearly the same number.
No mean to remove, and edge effects that shrink with scale. The Morlet wavelet has essentially no mean, and so nothing on this tab has to subtract a DC term or remove a per-window mean. The arbitrary horizontal origin simply never enters. And the untrustworthy region near the ends is no longer a straight band but a cone (shaded below) because a long coarse packet runs off the end of the data much sooner than a short fine one. That C-shaped opening is the wavelet’s version of the windowing effect from the previous tab, and it is worth reading carefully: the coarse scales along the left can only be measured well away from both ends, in the middle stretch of the tree, whereas the fine scales on the right stay usable almost to the trunk base and the branch tip.
What this does not fix. The cone eats any scale approaching the height of the tree you measured. That is a limit of the data and not of the method: a trajectory only a couple of bend-lengths tall cannot support a statement about bends, whatever transform you point at it. The heuristic divide is therefore useful only when it sits well below the height of the tree, which makes it a claim about how smooth the tree is rather than about where its fork happens to be.
x(z). Vertical axis is height, shared with the tree; horizontal axis is the
spatial frequency equivalent to each wavelet scale, on a log axis, as scales are
probed in constant ratios rather than constant steps. The shaded wedges are the
cone of influence, where the wavelet at that scale overruns the end of the data.
The dashed white line is the same heuristic divide fcut as the
previous tab.
One caution the smoothness of this picture hides: the transform is evaluated at every
height sample, and so the image is finer along the vertical axis than the measurement
behind it really is. The resolution of a row is the length of its packet, which makes it
different on every row, sharp along the fine scales at the right and very coarse along
the left, where neighboring columns rest on almost exactly the same stretch of
trajectory. The flare of the cone is that same length seen against the ends of the
record, and so the wider the cone is at a given row, the more the smoothness along that
row is telling you about the packet rather than about the tree.
hop. Note that this is
not a finer measurement, only a finer stride: neighboring packets overlap almost
entirely, and so the extra dots buy smoothness rather than information. What this tab
genuinely drops is the window length, not the slicing.
Morlet w0 — how many oscillations fit under the Gaussian envelope,
and the only shape parameter the wavelet has. Raise it and each packet carries more cycles,
which sharpens the range of frequencies it responds to at the cost of blurring where in
height it responds; lower it for the opposite trade. This is the same trade the previous
tab makes with W, with one difference that is the whole point of this tab:
here it applies to every scale at once rather than being set once for all of them. The
value 6 is the usual default, and it is what makes a packet’s scale and its wavelength
very nearly the same number, and so the scale axis can be read as a wavelength axis. Move it
to 3.8 and that correspondence stretches by about half, which is why the frequency range in
the readout shifts while the scale range does not.
Voices per octave — how many scales are computed per doubling of scale. More voices give a smoother scalogram at more compute and no more information, as neighboring scales overlap heavily, and so a blocky picture at a low setting is showing you the scale sampling rather than a shortage of pixels.
Lbend — the same heuristic length scale as the previous tab, and the
same divide fcut = 1/Lbend. It also decides which scales
feed the fine-scale profile on the right, namely everything to the right of the line.
Scale range — read-only, and worth reading. Neither end is a choice: the narrow end is two grid steps, the finest width the samples can represent, and the wide end is half the height of the tree, beyond which the cone covers every position and leaves nothing to measure. The count of scales then follows from those two ends and the voices setting.
log magnitude — rescale as log10(1+9v) so weak
components become visible. This is on by default here, unlike the previous tab,
because the tree’s slow bending dominates the coarse scales, and without it the finer
scales carrying the ant signal sit too far down the color scale to read.
per-slice normalization — rescale every height slice to its own maximum. This is off by default here, again unlike the previous tab, and it is worth leaving off. It destroys the comparison this tab exists to make, which is how much fine-scale energy one height carries relative to another. It is also worse than merely unhelpful: the deliberation bump lifts the coarse scales too, which raises each row’s maximum, and so the fine scales are pushed down the color scale exactly where the signal is. On the default trajectory it turns a contrast of 1.07 into 0.75, inverting it. The profile on the right is unaffected either way, being computed from unnormalized magnitudes.
shade cone of influence — show the region where the wavelet at that scale reaches past the end of the data. The cells inside it are not merely marked, as the edge bands on the previous tab are: they are left out of the picture, out of the profile, and out of the color normalization. Turn the shading off to see how much of the coarse end is resting on padding rather than on observations.
library(WaveletComp)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
# Uniform grid in HEIGHT, spaced to match the coarsest raw spacing in z so the
# scale range never claims resolution the samples do not have.
K <- round(diff(range(ants$z)) / max(diff(ants$z))) + 1
z_grid <- seq(min(ants$z), max(ants$z), length.out = K)
xz <- approx(ants$z, ants$x, xout = z_grid)$y
dz <- diff(z_grid)[1]
# analyze.wavelet works in steps of one sample, so periods come back in samples;
# multiply by dz for height.
#
# Three defaults of analyze.wavelet have to be overridden, and two of them were
# hiding each other. loess.span defaults to 0.75, which loess-detrends the
# series first; on this trajectory that removes a trend spanning 95% of the
# range of x, which is the tree shape this tab exists in order NOT to fit. The
# series is also zero-padded internally, and the ends here are nowhere near
# zero: the branch tip sits at x = 0.4, so the pad puts a step there. Setting
# loess.span = 0 alone therefore makes things look worse, not better, taking the
# top of the profile to 20x the interior. Hold the edge value out past the
# reach of the widest packet instead, then trim the padding back off.
# Finally, method defaults to "white.noise", whereas the null that matters here
# is red: an AR(1) trajectory has more coarse-scale power than white noise, and
# so "bright" is not the same as "surprising".
pad <- ceiling(4 * (K/2))
xpad <- c(rep(xz[1], pad), xz, rep(xz[K], pad))
wt <- analyze.wavelet(data.frame(x = xpad), "x",
loess.span = 0,
dt = 1, dj = 1/12,
lowerPeriod = 2, upperPeriod = K/2,
make.pval = TRUE, method = "AR",
params = list(AR = list(p = 1)), n.sim = 100)
keep <- (pad + 1):(pad + K) # drop the held-value margin again
Power <- wt$Power[, keep]
# Blank whatever the cone of influence reaches, as the widget does; a Morlet at
# scale s reaches sqrt(2)*s, in samples here as dt = 1. Those cells then draw as
# gaps instead of colour, and the profile below can average without them. The
# padding above is held values rather than data, so the cone is still needed;
# without it the top of the profile reads about 3x the interior.
reach <- sqrt(2) * wt$Scale
ok <- outer(reach, seq_len(K), function(r, i) (i - 1) >= r & (K - i) >= r)
Power[!ok] <- NA
# Plot the trimmed matrix directly. wt.image() would draw all 395 padded
# columns, of which only 79 are data, and it puts position on the horizontal
# axis; here VERTICAL = z, so it reads against a drawing of the tree.
freq <- 1 / (wt$Period * dz)
o <- order(freq)
image(x = freq[o], y = z_grid, z = Power[o, ],
log = "x", col = hcl.colors(64, "viridis"),
xlab = "spatial frequency (1/length, log)", ylab = "height z")
abline(h = 1, col = "darkorange", lwd = 2)
abline(v = 1 / 0.15, col = "white", lty = 2, lwd = 1.5)
# Fine-scale energy against height: the leading indicator. Everything shorter
# than L_bend counts as ant rather than tree.
L_bend <- 0.15
fine <- wt$Period * dz < L_bend
prof <- colMeans(Power[fine, , drop = FALSE], na.rm = TRUE)
plot(prof, z_grid, type = "l", xlab = "fine-scale power", ylab = "height z")
abline(h = 1, col = "darkorange", lwd = 2) # fork heightlibrary(WaveletComp)
library(dplyr)
library(tidyr)
library(tibble)
library(ggplot2)
# Inputs
# ants: one row per video frame, columns t (trajectory parameter running
# 0 to 1), x (lateral position), z (height up the tree).
K <- round(diff(range(ants$z)) / max(diff(ants$z))) + 1
z_grid <- seq(min(ants$z), max(ants$z), length.out = K)
xz <- approx(ants$z, ants$x, xout = z_grid)$y
dz <- diff(z_grid)[1]
# See the Base R tab for why loess.span, the padding and method are all set
# explicitly: the defaults detrend the series, zero-pad ends that are not near
# zero, and test against a white-noise null.
pad <- ceiling(4 * (K/2))
xpad <- c(rep(xz[1], pad), xz, rep(xz[K], pad))
wt <- analyze.wavelet(data.frame(x = xpad), "x",
loess.span = 0,
dt = 1, dj = 1/12,
lowerPeriod = 2, upperPeriod = K/2,
make.pval = TRUE, method = "AR",
params = list(AR = list(p = 1)), n.sim = 100)
# Trim the held-value margin, then blank whatever the cone of influence reaches
Power <- wt$Power[, (pad + 1):(pad + K)]
reach <- sqrt(2) * wt$Scale
ok <- outer(reach, seq_len(K), function(r, i) (i - 1) >= r & (K - i) >= r)
Power[!ok] <- NA
# Tidy the power matrix: one row per (height, scale)
scalo <- as_tibble(t(Power), .name_repair = "minimal") |>
setNames(as.character(wt$Period * dz)) |>
mutate(z = z_grid) |>
pivot_longer(-z, names_to = "scale", values_to = "power",
names_transform = list(scale = as.numeric)) |>
mutate(freq = 1 / scale)
ggplot(scalo, aes(freq, z, fill = power)) +
geom_raster() +
scale_x_log10() +
scale_fill_viridis_c(trans = "log10") +
geom_hline(yintercept = 1, color = "darkorange", linewidth = 1) +
geom_vline(xintercept = 1 / 0.15, color = "white", linetype = 2) +
labs(x = "spatial frequency (1/length, log)", y = "height z") +
theme_minimal()
# The same fine-scale profile, as a summary rather than a second plot.
# na.rm drops the cone cells blanked above; without it every height the cone
# touches comes back NA and the trace stops short at both ends.
scalo |>
# qualified because the previous tab's snippet attaches signal, whose filter()
# then masks dplyr's; a later library(dplyr) does not win the name back
dplyr::filter(scale < 0.15) |>
summarise(fine_power = mean(power, na.rm = TRUE), .by = z) |>
ggplot(aes(fine_power, z)) +
geom_path(color = "#6b21a8", linewidth = 1) +
geom_hline(yintercept = 1, color = "darkorange", linewidth = 1) +
labs(x = "fine-scale power", y = "height z") +
theme_minimal()