Labs on Queuing Theory, Signals, and Systems
Interactive labs I have built for my Cambridge course Computer Systems Modelling. Each one runs in the browser; no installation is needed. The cards below say what each lab shows, what to try, and which ideas from the course it illustrates.
Queuing theory
- Sampling by Inversion: drawing samples from any distribution
- Confidence Interval Lab: sample means, t-intervals and the CLT
- Random Walk: transience and recurrence in a birth–death chain
- Eigenvalue Lab: eigenvectors, stochastic matrices and the power method
- Markov Chain Lab: trajectories, distributions and stationarity
- CTMC Rate Lab: from coin flips to the rate matrix Q
- Birth–Death Lab: the continuous-time chains behind simple queues
- M/M/1 Lab: queue length, delay and Little’s theorem
- M/M/1/K Lab: finite buffers and blocking probability
- M/D/1 and M/G/1 Lab: service-time variability and the P-K formula
Signals and systems
- Phasor Lab: sinusoids as rotating vectors, and Euler’s formula
- Keying Lab: sending bits with ASK, FSK and PSK
- Convolution Lab: flip, slide, multiply and add
- Complex Exponential Lab: k est as a helix in three dimensions
- System Properties Lab: testing for linearity, time invariance and causality
- LTI Systems Lab: eigenfunctions, natural response and the impulse response
- Fourier Series Lab: periodic signals as sums of harmonics
- Fourier Transform Lab: transform pairs, properties and convolution
Queuing theory
Sampling by Inversion
Shows: how to generate samples from any distribution using the inverse-CDF method, X = F−1(U).
How it works: Sketch a density f with the mouse (or pick a preset: bimodal, piecewise constant, triangular, truncated exponential, uniform). The page computes the CDF F as you draw. Each sample picks u ~ U(0,1) on the vertical axis of F, runs across to the curve, and drops to x = F−1(u). A histogram of the samples builds up beneath, and converges to f.
Try this: Draw one sample at a time to follow the path, then switch off animation and draw batches. Note how steep parts of F catch many values of u while flat parts catch almost none.
Concepts: CDF and its inverse · inverse transform sampling · empirical distribution · sample mean vs. E[X]
Confidence Interval Lab
Shows: how the sample mean, sample variance and 95% confidence interval behave, and why the central limit theorem makes intervals work for any population.
How it works: In three steps. (1) Sketch a population density; its true mean μ and variance σ² are computed exactly. (2) Draw one run of n samples and watch X̄n, Sn² and the t-interval settle as samples arrive, with the interval narrowing like 1/√n. (3) Repeat for hundreds or thousands of runs: the sample means pile up as N(μ, σ²/n), and about 95% of the intervals contain μ.
Try this: Use a skewed population with small n (say 5) and see coverage fall short of 95%; then increase n and watch it recover. Compare a bimodal population with a uniform one: the histogram of means looks normal either way.
Concepts: sample mean and variance · Student’s t interval · coverage · central limit theorem · standard error
Random Walk
Shows: sample paths of a birth–death Markov chain: a random walk that steps up with probability p and down with probability 1 − p, held at its boundaries.
How it works: Choose p, a starting state, the number of steps, and either a finite chain on 1 … N or a semi-infinite chain on 1, 2, 3, …. Each press of Run adds a new trajectory in a new colour, so runs can be compared side by side. For the semi-infinite chain, a dashed line shows the expected drift s0 + (2p − 1)t.
Try this: On the semi-infinite chain, compare p = 0.4, 0.5 and 0.6. For p > ½ the walk escapes and tracks the drift line (transient); for p ≤ ½ it keeps returning to state 1 and the line is misleading (recurrent). On a finite chain, watch where the walk spends its time as p moves away from ½.
Concepts: Markov chains · reflecting boundaries · drift and linearity of expectation · transience vs. recurrence · √t spread of trajectories
Eigenvalue Lab
Shows: a refresher on eigenvalues and eigenvectors before they are used for Markov chains: what an eigenvector is, why every stochastic matrix has dominant eigenvalue 1, and how the power method finds the eigenvector that goes with it (§3.5–3.6).1
How it works: In three steps. (1) Drag a unit vector x around the circle and watch Ax: an eigenvector is a direction the matrix does not turn, and its eigenvalue is the stretch. The characteristic equation is solved alongside. (2) Edit a 3 × 3 stochastic matrix and check three facts: A1 = 1, Gerschgorin disks keep every eigenvalue inside the unit circle, and under mild assumptions 1 is the only eigenvalue on it. (3) Run the power method and watch the Rayleigh ratio converge to the dominant eigenvalue. Multiplying on the left gives the stationary distribution π.
Try this: Find both eigenvector directions of the book’s Example 3.15. Then try the rotation, which has none, and the deterministic 3-cycle, where three eigenvalues tie on the unit circle. Compare how fast the power method converges on Example 3.24 and on a chain where all states communicate. Each step has tasks with hidden answers.
Concepts: eigenvalues and eigenvectors · characteristic equation · stochastic matrices · spectral radius · Gerschgorin’s theorem · power method · Rayleigh ratio · left eigenvectors and the stationary distribution
Markov Chain Lab
Shows: a discrete-time Markov chain as a state graph and its stochastic matrix P, the difference between one random trajectory and the evolving distribution πt+1 = πtP, and two ways to find the stationary distribution.
How it works: In three steps. (1) Build a chain of 3–8 states from a preset or by editing P; the graph redraws as you type. (2) Step a single walker through the chain alongside the distribution πt, and compare both with 1,000 independent walkers and with the walker’s time average. (3) Find the stationary π as the dominant eigenvector of PT, by power iteration (with the Rayleigh quotient and the eigenvalues plotted in the complex plane) and by stepping through Gaussian elimination on (PT − I)πT = 0.
Try this: Load Two clusters, slow mixing and see how slowly power iteration converges when |λ2| is close to 1. Then try Period 2, where xk oscillates instead of converging, and Reducible, where the stationary distribution is not unique.
Concepts: stochastic matrices · trajectories vs. distributions · ensemble and time averages · stationary distribution · eigenvalues and power iteration · mixing rate · periodicity and reducibility
CTMC Rate Lab
Shows: how a continuous-time Markov chain arises from coin flips, what the entries of the rate matrix Q mean, and how Q relates to the transition probability matrix P of a discrete-time chain.
How it works: It builds up in steps. A Bernoulli trial, repeated until the first success, gives the geometric distribution, drawn as a two-state chain. Running one trial every Δt with p = qΔt and shrinking Δt turns the geometric into the exponential. With several states, a discrete chain that ticks every Δt runs beside a continuous one that jumps at any instant. Racing exponential alarm clocks show what qjk and qjj mean, and the page derives P(Δt) ≈ I + QΔt. Probability then flows between states like fluid between tanks, which leads to Equations 6.12 and 6.13.2
Try this: Drag Δt towards zero and watch the geometric staircase close in on e−qt. Run 2,000 clock races and check that the holding time is Exp(−qjj) and that each clock wins a share qjk/(−qjj) of the races. Then edit Q and press Jump to π* to see every tank’s inflow balance its outflow.
Concepts: Bernoulli, geometric and exponential distributions · memorylessness · rate matrix Q · P(Δt) = eQΔt · holding times · embedded chain · Equations 6.12 and 6.13 · stationary distribution
Birth–Death Lab
Shows: birth–death processes, the continuous-time Markov chains behind simple queues. From state j the process moves up at rate λj and down at rate μj (§6.4).3
How it works: In two tabs. (1) Finite chain: set λj and μj for up to ten states, or pick a preset (the book’s Example 6.18, M/M/1/K, an M/M/∞-style chain, or a bottleneck). The page draws the state-transition-rate diagram, simulates a sample path, and compares the time spent in each state with the stationary distribution π*j = π*0 ∏ λi/μi+1. (2) Semi-infinite chain: with the same λ and μ in every state (the M/M/1 queue), run several sample paths from one starting state, and watch the whole distribution π(t) evolve as a heat map, with E[N(t)] and π0(t) settling to ρ/(1 − ρ) and 1 − ρ.
Try this: On the finite chain, lengthen the horizon until the time fractions match π*. Set one λj to 0 and see the chain become reducible. On the semi-infinite chain, push ρ = λ/μ towards 1 and then past it: the paths stop forgetting where they started, E[N(t)] grows without bound, and the stationary distribution disappears.
Concepts: birth–death processes · state-transition-rate diagram · stationary distribution · time averages · M/M/1 and M/M/1/K queues · utilisation ρ · stability · transient behaviour of π(t)
M/M/1 Lab
Shows: the M/M/1 queue as a birth–death chain on 0, 1, 2, … with arrival rate λ and service rate μ in every state, its stationary distribution π*j = (1 − ρ)ρj, and the mean queue length and delay (§6.4.4).4
How it works: In four panels. (1) Sample paths from a chosen starting state, with the stationary mean ρ/(1 − ρ) as a reference line. (2) The distribution π(t), found by solving Equation 6.20 numerically, shown as a heat map and as a bar chart against π*, with E[N(t)] and π0(t) over time. (3) One long first-come-first-served run, whose running averages of customers and delay are compared with Equations 6.37 and 6.38 and checked against Little’s theorem. (4) A ρ slider that plots the M/M/1 formulas and the simulated averages together.
Try this: Start 15 customers deep and see how long the paths take to forget where they began. Try different initial distributions and watch π(t) reach the same limit. Move ρ from 0.5 to 0.95 and see delay grow like 1/(1 − ρ), and the long-run averages take longer to settle.
Concepts: M/M/1 queue · utilisation ρ · geometric stationary distribution · transient behaviour · mean number in system and in queue · mean delay · Little’s theorem · stability
M/M/1/K Lab
Shows: an M/M/1 queue with room for at most K customers (K − 1 waiting, one in service), where arrivals that find the system full are lost. Its stationary distribution π*j = (1 − ρ)ρj/(1 − ρK+1) gives the blocking probability π*K (§6.6.2, Equation 6.47).5
How it works: Set the utilisation ρ and the system size K, or load the book’s Example 6.25. Bars show π* beside the fraction of simulated time spent in each state. An animated queue shows buffers filling and arrivals being turned away, with running counts of arrivals and losses. A third panel tracks the fraction of arrivals lost and the fraction of time spent full: because Poisson arrivals see time averages, both converge to π*K.
Try this: Fix ρ = 0.8 and raise K to see how quickly the blocking probability falls. Set ρ above 1: the queue is still stable, but see where the probability sits. Run a few seeds and check whether the lost fraction stays above π*K, below it, or neither.
Concepts: finite buffers · M/M/1/K queue · blocking probability · throughput and loss rate · stability for ρ ≥ 1 · Poisson arrivals see time averages (PASTA)
M/D/1 and M/G/1 Lab
Shows: how the variability of service times affects delay in a single-server queue with Poisson arrivals, and the Pollaczek–Khinchine formula E(Wq) = (1/μ) · ρ(1 + CS²)/2(1 − ρ).
How it works: In two tabs. (1) M/D/1 against M/M/1: two animated queues share the same arrivals, one with exponential and one with constant service times. Their running mean delays are plotted against theory as you vary the service rate μ. (2) M/G/1 and the P-K formula: a second slider sets the squared coefficient of variation CS² of service time, from 0 (constant) to 4 (highly variable), using gamma-distributed service. The page works through the formula at the current settings and plots mean delay against both μ and CS². Either tab can show time in queue or time in system.
Try this: On the first tab, check that the M/D/1 queueing delay is exactly half the M/M/1 delay at every μ. On the second, snap between CS² = 0 and 1 to recover both cases, then push CS² to 4 and see how much longer the simulated mean takes to settle when rare, very long services dominate.
Concepts: M/D/1 and M/G/1 queues · service-time variability · squared coefficient of variation · Pollaczek–Khinchine formula · time in queue vs. time in system · utilisation ρ · convergence of simulated means
Signals and systems
Phasor Lab
Shows: how a sinusoid is the shadow of a rotating vector, and how this leads to the phasor c ejθ, the operator j as a quarter turn, Euler’s formula, and negative frequency (§5.1).6
How it works: In four tabs. (1) The rotating vector: a disk turns with period T; its projection on the X-axis traces out the carrier of Figure 5.1, and you can drag the vector by hand. (2) j as a rotation: multiply by j, by −1, or take the conjugate, and watch the vector turn or reflect. (3) Euler’s formula: edit a point in rectangular form a + jb or phasor form c ejθ and see the other follow, with buttons for Examples 5.1 and 5.2. (4) Negative frequency: two phasors turn in opposite directions, and their sum stays on the real axis (Equations 5.5 and 5.6). Each tab shows the matching Python code, with its output.
Try this: Reverse the direction of rotation and notice that the cosine does not change, but the sine does. Press Multiply by j four times and see where the vector ends up. Convert −3 + j4 to phasor form and work out why atan(b/a) gives the wrong angle. Each exercise has a hidden answer.
Concepts: sinusoids · phasors · the complex plane · j as rotation · Euler’s formula · rectangular and polar forms · conjugates · negative frequency
Keying Lab
Shows: how a transmitter sends bits by changing the phasor c ejωt from one symbol to the next: its magnitude (amplitude shift keying), its frequency (frequency shift keying), or its phase (phase shift keying), as in Example 5.3.7
How it works: Type a bit string and pick ASK, FSK, BPSK, QPSK, 8-PSK or 16-PSK. The page draws the transmitted waveform with the bits above each symbol, the rotating phasor being sent, and what the receiver sees after multiplying by its local oscillator: a constellation of points with decision boundaries, labelled with a Gray code. A noise slider adds random complex error to each symbol, and the decoded bits show errors underlined. A bar chart compares symbol error rates across the PSK schemes, and Python code for encoding, adding noise and decoding updates as you go.
Try this: Compare on-off keying with BPSK at the same noise and see why BPSK makes fewer errors. Switch to QPSK to send two bits per symbol. Raise the noise slowly and watch 16-PSK fail long before BPSK, which is why we cannot keep adding bits per symbol.
Concepts: carriers · amplitude, frequency and phase shift keying · BPSK and QPSK · constellations · decision regions · Gray codes · noise and symbol error rate
Convolution Lab
Shows: how to compute a convolution by flipping one function, sliding it along, multiplying point by point, and adding up the products, in both discrete and continuous time (§5.1.4–5.1.5).8
How it works: In two tabs. (1) Discrete time: type in x and y, or load the book’s Example 5.4, a delta, or a shifted delta. Five stacked plots show x, y, y flipped and slid to t, the products, and the output, with the sum written out term by term (Equation 5.11). (2) Continuous time: pick two functions, such as the triangle and rectangle of Example 5.5, and watch the overlap and its area trace out the output (Equation 5.16). A narrow unit-area pulse shows the Dirac delta in the limit. Python code for each case updates as you go.
Try this: Step t from 0 to 4 in Example 5.4 and check the values 3, 10, 25, 36, 69. Convolve with δ(t − 2) and see x shift. Load an asymmetric y and switch to Slide only to see why the flip matters. Shrink the pulse width and watch the output turn into a copy of x.
Concepts: discrete and continuous convolution · flip and slide · impulse function · Dirac delta · time shift · convolution vs. cross-correlation
Complex Exponential Lab
Shows: the complex exponential signal k est, with s = σ + jω, drawn in three dimensions, and how this one formula covers constants, real exponentials, sinusoids and damped sinusoids (§5.2.1).9
How it works: Set k, σ and ω with sliders, or by dragging s in the complex plane. The signal is drawn as a helix around the time axis, with its shadows on the Re–t and Im–t planes. Drag to turn the picture, or look straight down each axis. Presets reproduce the cases in Figures 5.7 and 5.8. A checkbox adds the conjugate k es*t, which makes the sum real. The two shadows are also shown flat, with their envelopes, and the page gives Python code to draw the same picture.
Try this: Start with σ = 0 and look down the t axis to see the helix turn. Make σ negative and then positive to see it shrink and grow. Change the sign of ω to reverse the direction of rotation. Set ω = 0 to flatten the helix into a real exponential, then add the conjugate to get a damped cosine.
Concepts: complex exponential · complex frequency s = σ + jω · the s-plane · growth and decay · damped sinusoids · conjugate pairs
System Properties Lab
Shows: how to test a system y(t) = H(x(t)) for linearity, time invariance, causality and memory by treating it as a black box, and why a failed test proves that a property is missing, while a passed test is only evidence (§5.3).10
How it works: In four tabs. (1) Linearity: pick a system, two inputs and two constants, and compare H(k1x1 + k2x2) with k1H(x1) + k2H(x2), with their difference plotted beneath. (2) Time invariance: compare delaying the input and then applying the system with applying the system and then delaying the output. (3) Causality and memory: add a small bump to the input and see where the output changes. (4) Classify them all: mark nine systems as causal, memoryless, time-invariant and linear, then check. Each tab shows the matching Python code, with its output.
Try this: Load Example 5.6 and see squaring fail superposition. Drive the saturation system with k1 = 0.5 and then 2 to see why limiting the output breaks homogeneity. Check that the offset y = x + 1 is not linear, even though its graph is a straight line. Compare the moving average of the last second with the centred average computed from a log.
Concepts: additivity and homogeneity · superposition · time invariance · causal and acausal systems · memoryless and dynamic systems · saturation and linearization · what a test can and cannot prove
LTI Systems Lab
Shows: the three questions of §5.4 for a linear time-invariant system: what it does to a complex exponential, what it does when the input is zero, and how to find its output for any input.11
How it works: In three tabs. (1) Exponential in, exponential out: feed k est into the systems of Examples 5.7 and 5.8 and watch input and output turn together in the complex plane. Their ratio stays constant, so the exponential is an eigenfunction (Equations 5.18–5.21). A plot shows the gain at each frequency. (2) Natural response: drag the roots of the characteristic polynomial, set y(0) and y′(0), and see the response of Equations 5.23 and 5.24. (3) Arbitrary input: cut the input into strips, add up a scaled, shifted impulse response for each, and watch the sum converge to the solution of the differential equation as the strips narrow (Equation 5.28).
Try this: Compare Example 5.7 with a squaring system to see why only an LTI system has an eigenvalue. Set s = j on Example 5.8 to find resonance. Move a root across the imaginary axis to turn decay into growth. Drive the lightly damped system at its natural frequency and watch the output build up, strip by strip.
Concepts: LTI systems · eigenfunctions and eigenvalues · transfer function · characteristic polynomial · natural frequencies and natural response · resonance · impulse response · convolution
Fourier Series Lab
Shows: how an eternal periodic signal is a sum of harmonics at multiples of the fundamental frequency, how each harmonic is a pair of rotating phasors, and why the coefficients of a pulse train lie on a sinc curve (§5.6).12
How it works: In three tabs. (1) Sum of harmonics: pick a square, triangle, sawtooth or pulse wave, or draw your own period, and add harmonics one at a time, with the coefficients of Equation 5.31 shown as a table and a bar chart. (2) Phasors and the complex form: see one harmonic as two phasors a quarter turn apart, drawn as helices in three dimensions, then a whole chain of phasors whose tip traces out the signal, in either the sinusoidal or the complex form (Equations 5.32 and 5.33). (3) Example 5.9: vary the width τ and period T0 of a pulse train and watch its coefficients move along the sinc envelope of Equation 5.36.
Try this: Add harmonics to the square wave and see that the overshoot at each jump stays near 9% however many you add (the Gibbs phenomenon). Compare how fast the square and triangle waves converge. Switch between the two forms and see that only the complex one keeps the tip on the real axis. Double T0 in Example 5.9 and see the samples crowd together while the zero crossings stay put.
Concepts: harmonics and the fundamental frequency · DC component · Fourier coefficients · even and odd signals · Gibbs phenomenon · complex form and two-sided spectrum · conjugate phasor pairs · the sinc function
Fourier Transform Lab
Shows: how the Fourier series of a pulse train becomes the Fourier transform of a single pulse as the period grows without limit, the standard pairs of Table 5.2, the properties that turn known transforms into new ones, and why convolution in time is multiplication in frequency (§5.7).13
How it works: In four tabs. (1) From series to transform: stretch the period of Example 5.10’s pulse train and watch the scaled coefficients fill in τ sinc(ωτ/2), with optional OFDM neighbours placed on its zeros. (2) Transform pairs: for each row of Table 5.2, damp the signal by e−ε|t| and watch its transform sharpen into impulses as ε → 0. (3) Properties: apply linearity, time shifting, time scaling, duality or differentiation and compare the predicted transform with one computed directly, in magnitude and phase. (4) Convolution: vary the system and input of Example 5.18 and compare the transform method with the convolution integral. Python code for each tab updates as you go.
Try this: Shrink ε for a constant and see 2πδ(ω) emerge. Slide a delay and see the magnitude stay put while the phase tilts. Halve the speed of a pulse and see its transform’s zero crossings move twice as close together. Check that differentiating a rectangular pulse gives two impulses whose transform is jω times the original.
Concepts: Fourier transform and its inverse · Dirichlet conditions · sinc and OFDM · impulses in frequency · linearity · time shifting and scaling · duality · differentiation · transfer function · the convolution property
Footnotes
Section numbers refer to Chapter 3 (Linear Algebra) of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Equation numbers refer to Chapter 6 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section, example and equation numbers refer to Chapter 6 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section and equation numbers refer to Chapter 6 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section, example and equation numbers refer to Chapter 6 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section, figure, example and equation numbers refer to Chapter 5 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Example numbers refer to Chapter 5 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section, example and equation numbers refer to Chapter 5 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section and figure numbers refer to Chapter 5 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section and example numbers refer to Chapter 5 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section, example and equation numbers refer to Chapter 5 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section, example and equation numbers refer to Chapter 5 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎
Section, example, table and equation numbers refer to Chapter 5 of S. Keshav, Mathematical Foundations of Computer Networking, Addison-Wesley, 2012.↩︎