Skip to content

Compressed Sensing: Recovering a Sparse Spectrum

Copy this iframe to your website:

1
<iframe src="https://dmccreary.github.io/signal-processing/sims/compressed-sensing/main.html" width="100%" height="708px" scrolling="no"></iframe>

Run MicroSim in Fullscreen

Description

The signal is \(N = 128\) samples long (1 s at 128 Hz) and real: \(x[n] = \sum_{i=1}^{k} a_i\cos(2\pi f_i n/N + \varphi_i)\), with \(k\) tones at whole-hertz frequencies from 1 to 63 Hz and amplitudes from 0.4 to 1. Because the signal is real, each tone occupies a conjugate pair of DFT bins, \(f\) and \(N - f\), so the spectrum \(\mathbf{x}_f\) is sparse. Only \(m\) samples are measured: \(\mathbf{y} = \Phi\mathbf{x}_f\) (plus noise), where \(\Phi\) holds the \(m\) rows of the inverse DFT matrix for the sample times \(t_i\).

Orthogonal matching pursuit (OMP) recovers \(\mathbf{x}_f\) in \(k\) iterations. Each iteration picks the frequency whose DFT column is most correlated with the residual, \(c(f) = \lvert\sum_i r_i e^{-j2\pi f t_i/N}\rvert\), adds both bins of its pair (a cosine and a sine column) to the support, and solves least squares on the support. The minimum-energy (least-squares) solution \(\Phi^H(\Phi\Phi^H)^{-1}\mathbf{y}\) is shown for contrast: it keeps the measured samples and sets the others to zero, which spreads its spectrum over many frequencies.

The top plot shows \(x(t)\) (blue), the samples (orange) and the OMP estimate (green); below are the minimum-energy spectrum and the OMP spectrum, with the true tones as blue rings. The dashed line at \(m/2\) Hz is the Nyquist frequency of \(m\) evenly spaced samples per second. Sliders set \(k\) (1 to 10), \(m\) (8 to 128) and the noise standard deviation \(\sigma\) (0 to 0.30). The menu chooses random or evenly spaced sample times, and New signal moves to the next seed (new tones, sample times and noise). Read-outs give both recovery errors \(\lVert\hat{x} - x\rVert/\lVert x\rVert\), whether the support was recovered, the \(m \ge C\,k\ln(N/k)\) rule, and the success rate over seeds 1 to 100. Reloading the page returns to seed 1.

This simulation is used in Chapter 14: Time-Frequency Analysis and Advanced Topics. It builds on the Sampling Theorem Explorer.

Lesson Plan

Learning Objectives

After using this MicroSim, students will be able to:

  1. Explain why \(m\) random samples can recover a \(k\)-tone spectrum whose tones lie far above \(m/2\) Hz, while \(m\) evenly spaced samples alias
  2. Analyze one run of OMP (correlate, select, least squares, update the residual) and relate its first pick to the tallest bar of the minimum-energy spectrum
  3. Compare the minimum-energy (least-squares) solution with the OMP solution using recovery error and sparsity
  4. Evaluate the rule \(m \ge C\,k\ln(N/k)\) against the success rates measured over 100 seeds, and estimate \(C\) for OMP
  5. Assess how measurement noise changes the recovery error and the support, and how extra measurements compensate

Target Audience

Undergraduate students in a junior or senior digital signal processing course. Plan 25 to 30 minutes.

Prerequisites

  • The DFT and its inverse, and the conjugate symmetry of the spectrum of a real signal
  • The sampling theorem and aliasing
  • Least squares for an overdetermined linear system

Activities

Activity 1: Recovery below the Nyquist rate (5 min). Use the defaults: seed 1, \(k = 3\), \(m = 24\) random samples (19%), \(\sigma = 0\). The tones are at 16, 42 and 56 Hz, all above \(m/2 = 12\) Hz. The banner is green, OMP picked 42, 16 and 56 Hz and the recovery error is 0.0%. The minimum-energy solution has a recovery error of 92.0%. Its tallest bar is the strongest tone (42 Hz), but bars at frequencies that are not tones (9, 58 and 28 Hz) are taller than the bars at 16 and 56 Hz. Discuss why OMP still finds 16 Hz second: once the 42 Hz tone is fitted and subtracted, the interference it spread over the other bins disappears from the residual.

Activity 2: How many samples? (8 min). Keep \(k = 3\) and lower \(m\) one step at a time. Seed 1 is recovered down to \(m = 15\) and fails from \(m = 14\) down. The 100-seed line gives the success rate:

\(m\) 12 16 20 24 28 32 36
support recovered (\(k = 3\)) 6 37 72 91 98 99 100

The rule of thumb predicts 23 to 45 samples. Repeat with \(k = 5\) (rule: 32 to 65): 38, 65, 82, 95, 99 and 100 of 100 at \(m\) = 24, 28, 32, 36, 40 and 48. The success rate first reaches 95% at \(m = 28\) for \(k = 3\) and at \(m = 36\) for \(k = 5\); the read-out shows these as \(C = 2.5\) and \(C = 2.2\), inside the rule's range.

Activity 3: Random or evenly spaced? (7 min). Reload the page and set \(m = 32\). With random times, seed 1 is recovered and the 100-seed line reads 99 of 100. Choose Evenly spaced times: now every 4th sample is taken, 32 samples/s with a Nyquist frequency of 16 Hz. The minimum-energy spectrum shows equal copies of each tone (the 42 Hz tone at 10, 22, 42 and 54 Hz; the 56 Hz tone at 8, 24, 40 and 56 Hz). OMP picks 10, 16 and 8 Hz: the aliases \(|42 - 32| = 10\) and \(|56 - 64| = 8\), exactly as in the Sampling Theorem Explorer. The 16 Hz tone sits exactly at \(m/2\), so its stem stops just below its ring (0.45 instead of 0.47). Only 2 of 100 seeds are recovered. At \(m = 64\) evenly spaced, OMP picks 22, 16 and 8 Hz, and 12 of 100 seeds are recovered.

Activity 4: Noise (5 min). Reload, set \(m = 32\) and raise \(\sigma\). For seed 1 the recovery error is 5.9% at \(\sigma = 0.10\) and 11.7% at \(\sigma = 0.20\), and the support is still recovered. At \(\sigma = 0.30\) OMP picks 31 Hz instead of 56 Hz and the error jumps to 50.7%. Over 100 seeds: 99 of 100 (median error 5.0%) at \(\sigma = 0.10\), 99 (9.9%) at 0.20 and 93 (14.9%) at 0.30. Now raise \(m\) to 48 with \(\sigma = 0.30\): the support is recovered, the error is 12.9%, and all 100 seeds succeed.

Assessment

  1. For \(k = 5\) tones and \(N = 128\), what range of \(m\) does \(m \ge C\,k\ln(N/k)\) with \(C\) = 2 to 4 predict, and how does the sim compare? Answer: \(k\ln(N/k) = 5\ln 25.6 = 16.2\), so 32 to 65 samples. The sim recovers 82 of 100 seeds at \(m = 32\), 95 at \(m = 36\) and all 100 at \(m = 48\).
  2. Why is the error of the minimum-energy solution close to \(\sqrt{1 - m/N}\)? Answer: by Parseval's theorem, minimizing the spectrum's energy minimizes the signal's energy, so the unmeasured samples are set to zero. The error is the energy of the missing samples, on average a fraction \(1 - m/N\) of the total. Seed 1 with \(k = 3\): 88.6% at \(m = 32\) (\(\sqrt{0.75} = 86.6\%\)) and 69.6% at \(m = 64\) (\(\sqrt{0.5} = 70.7\%\)).
  3. With 32 evenly spaced samples per second, OMP reports a 42 Hz tone as 10 Hz. Why do 32 random sample times not have this problem? Answer: on the uniform grid the DFT columns for \(f\), \(32 - f\), \(32 + f\) and \(64 - f\) have identical correlations with any real residual, so the 42 Hz tone is indistinguishable from 10 Hz. Random times make every column different; the aliasing becomes low-level interference spread over all frequencies, and the true tones stand out.
  4. Can \(m = 10\) samples recover \(k = 5\) tones? Answer: no. Five tones have 15 unknowns (5 frequencies, 5 amplitudes and 5 phases), more than the 10 equations; the sim recovers 0 of 100 seeds. In practice \(m \approx 2.2\,k\ln(N/k) = 36\) samples give 95% success.
  5. With the defaults, why is OMP's first pick 42 Hz? Answer: the first correlation is \(c(f) = \lvert\sum_i y_i e^{-j2\pi f t_i/N}\rvert\), which is \(N/2\) times the height of the minimum-energy bar at \(f\) (for 1 to 63 Hz). The tallest bar is at 42 Hz, the strongest tone (\(a = 0.92\)).

Common Misconceptions

"Compressed sensing beats the sampling theorem for any signal." It needs sparsity. With \(k = 10\) and \(m = 32\) only 1 of 100 seeds is recovered; with \(m = 64\) all 100 are. A signal with many tones needs many samples.

"Fewer samples just means a lower uniform sampling rate." Evenly spaced samples alias: at \(m = 32\) OMP returns the alias frequencies 10 and 8 Hz and only 2 of 100 seeds succeed. The same number of random samples succeeds for 99 of 100 seeds.

References

  1. Compressed sensing - Wikipedia - sparsity, incoherent measurements and the \(m \ge C\,k\log(N/k)\) result
  2. Matching pursuit - Wikipedia - the greedy algorithm and its orthogonal variant (OMP) used in this sim
  3. Restricted isometry property - Wikipedia - the condition on \(\Phi\) that guarantees sparse recovery
  4. Nyquist–Shannon sampling theorem - Wikipedia - the uniform-sampling limit that random sampling of a sparse spectrum gets around
  5. Sparse approximation - Wikipedia - finding a few columns of a dictionary that explain a measurement
  6. The Scientist and Engineer's Guide to Digital Signal Processing - Steven W. Smith - free online textbook on the DFT, sampling and aliasing