CASE STUDY · 03 · ENGINEERING COMPUTATION

SIMULATION, SIGNALS & LEARNING

Three computation-heavy projects, one habit: build the mathematics yourself. A shaker-table test rig driven through a prescribed 3D dance by three cam-driven followers; the signal-processing tools built along the way — a Kalman filter, displacement from an accelerometer, a Fourier audio compressor; and an early machine-learning study reading Alzheimer's off a brain signal.

Z = Z + X sin(ψ) cos(θ) Y sin(θ)  ·  ψ = 0.25 e⁻¹⁰ᵗ cos²t  ·  θ = 0.5 e⁻¹⁰ᵗ sin²t
SCROLL TO ENTER
01  THE PROJECT · DYNAMICS OF MACHINERY

A platform with three programmed degrees of freedom

This was my Dynamics of Machinery course project: design a shaker test table — a rigid platform driven along a prescribed trajectory to qualify a payload against shock and fatigue.

The platform rests on three followers. Each follower is lifted by a cam through a ball-screw, so three independent lift schedules give the platform three degrees of freedom: roll about Y (ψ), pitch about X (θ), and heave along Z. The target trajectory is a damped transient whose amplitudes are fixed by my student number — 5040, mapped digit-by-digit to (a, b, c, w) with 0 → 10, giving (5, 10, 4, 10) and λ = −w = −10:

ψ(t) = 0.25 e⁻¹⁰ᵗ cos²(t)   // roll, about Y
θ(t) = 0.50 e⁻¹⁰ᵗ sin²(t)   // pitch, about X
Z(t) = 0.20 e⁻¹⁰ᵗ sin²(t)   // heave, along Z

Every later result — the cam profiles, the joint accelerations, the cam and jack forces — is derived from these three laws. The kinematics were built symbolically in SymPy and SciPy, and the assembled mechanism was verified in MSC Adams. The panel evaluates the exact pose at any instant t ∈ [0, 0.5] s.

SHAKER TABLESYMPYMSC ADAMSJune 2023
LIVE · PLATFORM POSE ψ, θ, Z
t (s)
ψ · roll
θ · pitch
Z · heave
Rectangle 2 m × 1 m on three followers. Corners transformed by Ry(ψ)·Rx(θ) then heaved by Z; heave exaggerated ~2× for legibility.
READ MORE · MECHANISM & PARAMETERSCAM → FOLLOWER → BALL-SCREW

The rig converts three rotating cams into a fully prescribed rigid-body motion of the platform. Each actuation stage is deterministic, so the platform pose is a closed-form function of time — the property that makes the whole analysis analytic rather than numerical.

01 · CAMRotating profile sets the lift schedule
02 · FOLLOWERWheel, r = 5 cm, rides the profile
03 · BALL-SCREWConverts lift to vertical joint travel
04 · PLATFORMThree lifts → ψ, θ, Z

Why three followers give three DOF

Three non-collinear support points fully constrain a rigid plane. Prescribing the vertical position of each support therefore prescribes the platform's orientation and height uniquely: the two independent tilts (roll ψ, pitch θ) and the common heave Z. The three cams are the control inputs; the platform pose is the output.

The parameters

Amplitudes

(a, b, c) = (5, 10, 4) scale ψ, θ, Z. The digit map (student ID, 0 → 10) is the assignment's way of giving every student a distinct trajectory.

Decay & frequency

w = 10 sets the oscillation rate (0.1·w = 1 rad/s) and λ = −w = −10 the decay, so the transient dies within ≈ 0.5 s.

Physical constants

Wheel radius 0.05 m, jack speed 1 m/s, platform mass 10 kg, analysis window t ∈ [0, 0.5] s.

Because every input is a smooth, decaying harmonic (note the sin²/cos² envelopes, which start and end with zero slope), the platform executes a single controlled shudder and returns to rest — exactly the transient a shock-qualification test applies.

02  FINDING THE POSITION FUNCTION

The position function

Every downstream quantity depends on a single expression: the vertical position of each follower foot as a function of time. I derive it with 4×4 homogeneous transformation matrices, the standard construction for composing rigid-body motions.

A point is carried as the row vector [X Y Z 1]. In homogeneous coordinates a translation becomes a matrix multiply, so a rotation about Y, a rotation about X, and a translation compose into one transition matrix by ordinary multiplication, and the transformed position is read directly off the product:

Transition = Ry(ψ) · Rx(θ) · T
[X  Y  Z  1] · Transition = [X  Y  Z  1]

With the in-plane translations set to zero and a point taken on the platform surface (Z0 = 0), the vertical coordinate reduces to one closed form:

Z(t) = Z(t) + X·sin(ψ)·cos(θ) Y·sin(θ)

Substituting ψ(t), θ(t), Z(t) makes Z a pure function of time. Evaluated at the three feet — joint 1 (−1, −0.5), joint 2 (0, 0.5), joint 3 (1, −0.5) — it is the exact lift each cam must produce. Joint 2 lies on the roll axis and barely moves; joints 1 and 3 move hardest and in opposition. The complete matrix derivation, exactly as solved, is expanded below.

HOMOGENEOUS TRANSFORMSKINEMATICSSYMPY
LIVE · Z-DISPLACEMENT OF EACH FOLLOWER
t (s)
z₁
z₂
z₃
The exact curves the SymPy model produced — a damped transient that dies away in half a second.
READ MORE · THE FULL DERIVATIONRy·Rx·T, EXPANDED

This is the derivation exactly as solved in SymPy — the three base transforms, their product, the transformed position, and the reduction to a single scalar height function. A point is written as a row vector and multiplied on the left, so the composed transform reads left-to-right.

1 · The three base transforms

A rotation about Y by ψ, a rotation about X by θ, and a pure translation, each a 4×4 matrix acting on [X Y Z 1]:

Ry(ψ) =cos ψ0sin ψ00100−sin ψ0cos ψ00001
Rx(θ) =10000cos θ−sin θ00sin θcos θ00001
T(X, Y, Z) =100001000010XYZ1

2 · The transition matrix

Multiplying Transition = Ry · Rx · T gives the single 4×4 that maps a rest position to its pose at time t:

Transition =cos ψsin ψ sin θsin ψ cos θ00cos θ−sin θ0−sin ψsin θ cos ψcos ψ cos θ0XYZ1

3 · The transformed position

With [Xpos Ypos Zpos 1] · Transition = [X′ Y′ Z′ 1], the three transformed coordinates are:

X′ = X + Xpos cos ψ − Zpos sin ψ
Y′ = Y + Xpos sin ψ sin θ + Ypos cos θ + Zpos sin θ cos ψ
Z′ = Z + Xpos sin ψ cos θ − Ypos sin θ + Zpos cos ψ cos θ

4 · Reduction for this rig

There is no in-plane translation (X = Y = 0) and each foot is a point on the platform surface (Zpos = 0). The in-plane outputs then satisfy X′ ≈ Xpos, Y′ ≈ Ypos by the geometry of the design, and the height collapses to:

Z′pos = Z + Xpos sin ψ cos θ − Ypos sin θ

5 · Substituting the trajectory (a, b, c, w) = (5, 10, 4, 10)

ψ = 0.25 e^(−10t) cos²(t)  ·  θ = 0.5 e^(−10t) sin²(t)  ·  Z = 0.2 e^(−10t) sin²(t)

6 · The height function, in full

Combining steps 4 and 5 gives the exact scalar the whole project runs on — height as a function of a foot's rest position (Xpos, Ypos) and time:

Z′pos = Xpos sin(0.25 e^(−10t) cos²t) · cos(0.5 e^(−10t) sin²t) − Ypos sin(0.5 e^(−10t) sin²t) + 0.2 e^(−10t) sin²t

Why this matters: because Z′ is closed-form, velocity and acceleration are exact derivatives, the cam profiles are exact envelopes, and the contact forces are exact — nothing downstream is a numerical approximation.

03  DESIGNING THE CAM PROFILES

Cam profiles as the envelope of the follower

The follower contacts the cam through a wheel of radius 5 cm, not at a point. The machinable cam surface is therefore not the lift curve itself but the envelope of the family of wheel circles swept along it.

For each follower I take its lift z(s) — where s is the traversed cam coordinate (jack speed 1 m/s, so s = t) — place a 5 cm circle at every sample, and take the boundary of their union. That envelope is the profile that will be cut. The three profiles differ sharply: cam 1 falls and hooks, cam 2 is nearly straight because joint 2 hardly moves, and cam 3 is the mirror of cam 1. Select a cam and scrub the wheel to see it ride the profile it generates.

ENVELOPE OF CIRCLESCAM PROFILEr = 5 cm
LIVE · CAM PROFILE FROM THE ENVELOPE
CAM
WHEEL
follower1
lift z
radius5 cm
Black body = union of 5 cm wheel circles along the track. Dashed line = the follower's displacement (the wheel-centre path).
READ MORE · THE ENVELOPE CONSTRUCTIONFROM LIFT CURVE TO CUT PROFILE

The wheel centre follows the lift curve, but the surface the cam is cut to is offset from it by the wheel radius. The correct profile is the envelope of the one-parameter family of wheel circles centred on that curve.

The envelope condition

Let the wheel centre trace (s, z(s)). The family of contact circles is

F(x, y, s) = (x − s)² + (y − z(s))² − r² = 0

The envelope is the curve tangent to every member of the family, found by solving F = 0 together with ∂F/∂s = 0 and eliminating s. That yields the two parallel (offset) curves at normal distance r from the lift curve:

xenv = s ∓ r·z′(s) / √(1 + z′(s)²)  ,  yenv = z(s) ± r / √(1 + z′(s)²)

The lower offset is the working face the follower rides. Where the lift is flat (z′ = 0) the offset is simply r below the curve; where it steepens, contact shifts along the wheel and the offset rotates with the local normal — which is why the cam face bulges away from the dashed centre-line on the curved sections.

How the panel computes it

Discrete union

Rather than solve the elimination symbolically, the profile is rendered as the union of the wheel circles at every sampled s — the same numerical envelope used in the report. The outer boundary of that black region is the cut profile.

Why per-cam shapes differ

The envelope inherits the lift curve's shape: joint 2's near-zero lift gives an almost straight bar; joints 1 and 3 produce opposite hooks because their lifts are mirror transients.

04  MASS DISTRIBUTION, NORMAL & FORCES

Contact forces from the free body

With the motion in closed form, the loads follow from Newton's second law resolved along the cam normal. The supported mass is fixed first: enforcing the total mass and the two centre-of-mass conditions distributes the 10 kg platform as 2.5 / 5 / 2.5 kg over the three joints.

Differentiating Z twice gives each follower's vertical acceleration az. The contact force acts along the cam surface normal, whose inclination is θn = atan2(1, −vz). Resolving the free body vertically and horizontally yields the force on the cam and the force each jack must supply to hold constant speed:

Fcam = (mj + m)(g + az) / sin θn
Fjack = (mj + m)(g + az) · cot θn

The demand peaks at t = 0 — joint 1 requires ≈ 260 N in the first instant — then decays with the transient, while joint 2 stays near zero. Toggle jacks and cams to compare; the full free-body derivation is expanded below.

FREE-BODYa_z = d²Z′/dt²CONTACT FORCE
LIVE · REQUIRED FORCES vs TIME
SHOWt
F₁ (N)
F₂ (N)
F₃ (N)
Reproduces the report's force diagrams: ≈ 3.5 kg supported per arm, g = 9.81, forces resolved through the surface normal.
READ MORE · MASS BALANCE & FREE-BODYSTATICS + NEWTON ALONG THE NORMAL

The loads split into two problems: a static one that fixes how much mass each joint carries, and a dynamic one that resolves the inertial and gravitational load through the inclined cam contact.

1 · Mass distribution

Three unknowns m1, m2, m3; three equations — the total mass and the two centre-of-mass conditions that keep the platform balanced over its supports at (Xpos, Ypos):

m1 + m2 + m3 = 10
X1m1 + X2m2 + X3m3 = 0  ,  Y1m1 + Y2m2 + Y3m3 = 0

With the feet at (−1, −0.5), (0, 0.5), (1, −0.5), the solution is symmetric:

(m1, m2, m3) = (2.5, 5, 2.5) kg

2 · The surface normal

The follower rides the cam face; the reaction acts along its normal. From the profile slope vz = dZ′/dt, the normal inclination is

θn = atan2(1, −vz) = arctan( 1 / −vz )

3 · Free-body equations

For an arm carrying total mass mj + m ≈ 3.5 kg (the arm plus its share of the platform), moving at acceleration az, vertical and horizontal balance give:

ΣFz:  (mj+m) az = −(mj+m) g + Fcam sin θn  ⇒  Fcam = (mj+m)(g + az) / sin θn
ΣFx:  Fjack = Fcam cos θn = (mj+m)(g + az) · cot θn

The jack only has to overcome the horizontal component of the cam reaction to keep the follower advancing at constant speed, so it scales with cot θn; the cam itself carries the full inclined reaction, scaling with 1/sin θn. Both are evaluated live from the exact az(t) and vz(t), reproducing the report's force diagrams.

05  CONTROL · STATE ESTIMATION

A Kalman filter, implemented from scratch

This one came out of my Control course, which I studied and implemented on my own. A Kalman filter is the optimal linear estimator: it fuses a noisy measurement stream with a model of the system's motion to recover the state.

I implemented the constant-velocity form. The state carries position and velocity; each step predicts forward with the motion model, then corrects against the new measurement, weighting the two by their covariances through the Kalman gain. The shaded band is the filter's own ±2σ confidence, which contracts as evidence accumulates. Raising the measurement noise R makes the estimate lean on the model; raising the process noise Q makes it track the data more aggressively. The full recursion is below.

KALMAN FILTERPREDICT · UPDATEOPTIMAL ESTIMATOR
LIVE · TRUTH · MEASUREMENTS · ESTIMATE
MEAS NOISE R
PROCESS NOISE Q
raw RMSE
filtered RMSE
gain K
A 2-state (position, velocity) linear Kalman filter running live in the browser.
READ MORE · THE KALMAN RECURSIONPREDICT / UPDATE, EXACTLY

For a linear system with Gaussian noise, the Kalman filter is provably the minimum-variance estimator. It maintains a state estimate and its covariance, and updates both every step in two stages.

Model

Constant-velocity motion: the state is position and velocity, propagated by F, observed in position by H.

x = [ p, v ]ᵀ  ,  F = [ [1, dt], [0, 1] ]  ,  H = [ 1, 0 ]

Q is the process-noise covariance (how much the model is trusted); R is the measurement-noise variance.

Predict

x⁻ = F x  ,  P⁻ = F P Fᵀ + Q

Update

y = z − H x⁻  (innovation)
S = H P⁻ Hᵀ + R  ,  K = P⁻ Hᵀ S⁻¹  (Kalman gain)
x = x⁻ + K y  ,  P = (I − K H) P⁻

What the gain means

The R / Q ratio sets the trust

K is large when the prediction is uncertain relative to the measurement (small R, large Q) and small in the reverse case. The estimate is a running, covariance-weighted average of model and data.

Implementation detail

The filter is seeded with velocity from the first two samples and P initialised to the measurement variance, so it converges without a large startup transient — the tuning that fixed the earlier over-noisy behaviour.

06  SIGNAL PROCESSING · DOUBLE INTEGRATION

Displacement from an accelerometer

Recovering displacement from an accelerometer is a double integration: velocity is the integral of acceleration, displacement the integral of velocity. Numerically, each integral is the running area under the curve.

The panel makes that concrete — a sweeping cursor accumulates the signed area under the acceleration curve into the velocity, then the area under velocity into displacement, with the running integral shown at each step.

The difficulty is bias. A constant sensor offset integrates to a linear velocity ramp and a quadratic displacement drift, so the naive result diverges. The remedy is spectral, not mechanical: remove the low-frequency component — a detrend / high-pass — between integrations. Toggle drift correction to see the recovered displacement return to the true motion. The mechanism is expanded below.

∫ a dt = v∫ v dt = xHIGH-PASS · DETREND
LIVE · a → v → displacement
SENSOR BIAS
NOISE
peak drift
disp error
state
Top: measured acceleration (with bias). Middle: velocity. Bottom: recovered vs true displacement.
READ MORE · WHY BIAS BECOMES DRIFTA CONSTANT OFFSET → A PARABOLA

The whole problem is that integration accumulates error. A tiny bias that is invisible in the raw acceleration becomes the dominant term after two integrations.

The error propagation

Suppose the measurement is the true acceleration plus a constant bias b:

ameas(t) = a(t) + b
v(t) = ∫ ameas dt = ∫ a dt + b·t  (a linear ramp)
x(t) = ∫ v dt = ∫∫ a dt + ½ b·t²  (a quadratic drift)

So a constant offset produces a parabola that swamps the real displacement — exactly the divergence you see with correction off.

The fix, and why it works

Separation in frequency

The drift lives at very low frequency (DC and near-DC); the vibration of interest sits above it. A high-pass / detrend step removes the low band without touching the signal band.

Implementation

Integration is trapezoidal; after each stage the series is de-meaned and a least-squares linear trend is subtracted, cancelling the bias-induced ramp before it can compound into the next integral.

This is the standard pipeline for reading a vibrating structure's own accelerometers — the sensors a shaker table like the one in this case study actually carries.

07  ENGINEERING MATHEMATICS · FOURIER COMPRESSION

Audio compression in the frequency domain

This was my Engineering Mathematics project. A signal is a sum of frequencies, and in the frequency domain most coefficients carry negligible energy — so discarding them compresses the file with little audible loss.

The compressor takes a short-time Fourier transform, keeps only the highest-energy coefficients per frame, and reconstructs from those alone. Retaining the largest coefficients is the optimal choice for a fixed budget — the best L2 approximation of the frame — so the result is as close to the original as any method could make it at that size. As the budget drops the spectrum thins, the file shrinks in real time, and the SNR quantifies the loss.

Press play to compare the original and the compressed reconstruction directly. The derivation and the size model are expanded below.

FOURIER · STFTTOP-K COEFFICIENTSA / B LISTEN
LIVE · FOURIER COMPRESSION · HEAR THE DIFFERENCE
KEEP COEFFS
size
saved
SNR
Top: original (muted) vs Fourier-reconstructed (accent) wave, with a playback cursor. Bottom: one frame's spectrum — kept coefficients in accent, discarded in grey, threshold dashed.
READ MORE · STFT, OPTIMALITY & SIZEWHY TOP-K IS THE BEST YOU CAN DO

The method is transform coding — the principle behind MP3 and JPEG. Move to a basis where energy concentrates, keep the few coefficients that matter, and store only those.

01 · FRAMEHann window, 50% overlap
02 · FFTPer-frame spectrum
03 · KEEP TOP-KLargest coefficients
04 · INVERSEIFFT + overlap-add

Why the top-K coefficients are optimal

The Fourier basis is orthonormal, so by Parseval's theorem the squared reconstruction error equals the sum of the squared magnitudes of the discarded coefficients:

‖ x − x̂ ‖² = Σdiscarded |ck

To minimise that error while keeping exactly K coefficients, you must discard the K-smallest in magnitude — i.e. keep the K-largest. No other K-term choice in this basis can do better; that is the sense in which the compression is optimal.

The size model

Raw audio costs samples × bit-depth. A kept coefficient must store its index and its complex value, so for a frame of length F keeping a fraction p costs:

saved ≈ 1 − p · ( log₂F + 2·bval ) / (bit-depth)

With F = 1024 and 12-bit values against 16-bit PCM, keeping ~10% of coefficients saves ~85–90%. Fidelity is reported as the signal-to-noise ratio over the whole clip:

SNR = 10 · log₁₀ ( Σ x² / Σ (x − x̂)² )  [dB]

Note: the reconstruction uses overlapping windows for audio quality, while the size figure is reported against a critically-sampled coefficient count — the standard accounting for a transform codec.

08  WHERE IT STARTED · HIGH SCHOOL · 2018

Reading Alzheimer's off a brain signal

My first machine-learning project, written in 2018 when I was fifteen and still in high school. The question: can a model detect Alzheimer's disease earlier and more reliably than routine clinical diagnosis?

All of the data and figures here are from that 2018 study. Clinical diagnosis then sat near 87% sensitivity and 58% specificity (median). Working from event-related potential (ERP) markers extracted from EEG, the pipeline applied pre-processing and feature extraction, a recursive feature-selection loop, and an electrode-weighting step to classify subjects as AD or non-AD. It reached 85–91% accuracy, above the prior machine-learning benchmark of ≈ 86%. It was my entry point into machine learning, and into turning noisy real-world signals into a decision.

EEG · ERP MARKERSFEATURE SELECTION2018
RESULT · DIAGNOSTIC ACCURACY · 2018
Median clinical accuracy vs. the prior ML best vs. this study's ERP-marker classifier (85–91%). Figures from the 2018 study.
READ MORE · THE 2018 STUDYERP PIPELINE & RESULT

All figures and references in this section are from the original 2018 work. Alzheimer's is a slowly progressing neurodegenerative disease; earlier, more reliable detection was — and remains — clinically valuable, and the study asked whether EEG-derived ERP markers plus machine learning could improve on routine diagnosis.

01 · ERP / EEGEvent-related potential markers
02 · FEATURESPre-processing + extraction
03 · SELECTIONRecursive feature selection
04 · CLASSIFYElectrode weighting → AD / non-AD

Why ERP markers

Raw EEG is high-dimensional and hard to use clinically. Event-related potentials give a cleaner, time-locked response of the brain to a stimulus, which makes the signal tractable for a classifier while still carrying the disease-relevant information.

Result

Baseline (2018)

Median clinical diagnosis: 87% sensitivity, 58% specificity. Prior ML best: ≈ 86%.

This study

ERP-marker classifier: 85–91% accuracy, clearing the target set at the outset.

Context

A high-school research project (2018) — my first, and the reason the later signal-processing and simulation work felt natural.

·  SUMMARY

Four projects, one method

A cam-driven shaker table, a Kalman estimator, a Fourier audio codec, and an ERP classifier span mechanical dynamics, control, signal processing, and machine learning. Each was derived from its governing mathematics and implemented directly, so every result on this page is reproducible from the equations it rests on.