In [2]:
import numpy
import CS
from numpy import pi
from numpy.random import randn, random_sample
from IPython.html import widgets
from IPython.display import display
from matplotlib import pyplot as plt
# Initialize the demo
demo = CS.demo()
%matplotlib inline
Compressed sensing (CS) is an emerging field of mathematics and engineering that challenges the conventional paradigms of data acquisition. Since the first publication(Candes, 2007), the field has developed a substantial academic literature and has provided the foundation for major innovations in medical imaging, astronomy, and digital photography. Although the field of CS is very mathematical, it is conceptually intuitive. This tutorial provides a brief introduction to CS through simple examples and interactive figures.
Many signals have an inherent low-dimensional structure, meaning all the information can be compressed into relatively few coefficients. As a simple example, consider the signal $x(t)=Asin(ft +\phi)$ which is constructed with only three pieces of information: the amplitude $A$, the frequency $f$, and the phase $\phi$. Measuring this signal using conventional acquisition requires sampling the entire duration of the signal at a rate of at least $2f$, which for a 10 second 30 Hz signal results in 600 samples. Why did we collect 600 samples of data when we only require 3 pieces of information? This is the fundamental question CS is asking, and in a field with large data acquisition costs we should pay close attention to the answers.
This tutorial demonstrates the basic concepts of CS to the acquisition of the simple 1D signal:
$x(t) = 50\sin(2\pi f_1 t + \phi_1) + 100\sin(2\pi f_2t + \phi_2)$
In [1]:
xt, t = demo.test_signal()
The first requirement for compressed sensing is the existence of a sparse signal representation. We need to know a priori that the signal we are acquiring has relatively few non-zero coefficients in some transform domain. This may seem like a heavy requirement, but fortunately the field of data compression has found sparsifying transforms for many types of general signal classes. For example, the wavelet transform sparsifies images, the Fourier transform sparsifies harmonic signals, and curvelets sparsify seismic data. Data compression formats such as mp3 and jpeg rely on these sparse representations to reduce the size of files by keeping only the largest coefficients. For our simple example, let's assume we know the signal has harmonic content so it will be sparse in the Fourier domain.
In [4]:
plt.subplot(211)
demo.ts_plot(t, xt)
plt.subplot(212)
Xf,f = demo.fft(xt, t)
demo.fourier_plot(f,Xf)
Uniform sampling
Uniform sampling has been the engineered data acquisition strategy for decades. Shannon-Nyquist theory says we can uniquely represent our signal by sampling at least twice the highest frequency, which can be an expensive requirement. Considering our test signal can be reduced to only 2 samples of information, we are grossly oversampling. Let's see what happens when try to collect less using uniform sampling.
In [5]:
sample_slider = widgets.IntSliderWidget(min=1, max=9, step=1, value=6)
w=widgets.interactive(demo.uniform_sample,downsample_factor=sample_slider)
sample_slider.msg_throttle = 1
display(w)
Notice we can no longer reconstruct our signal as we decrease the number of samples below the Nyquist rate. The extra frequency coefficients that creep into spectra are called aliases, and they set the lower limit on the number of samples we need to acquire in order to uniquely define our signal.
Aliasing
Let's look closer at the cause of aliasing. Fundamentally, aliasing occurs when multiple signals have the same value at every sample, which creates ambiguity to what the true signal is. The figure below demonstrates aliased signals; the two signals appear identical once they are sampled.
In [6]:
demo.aliased()
The problem is that both the signal and sampling pattern are periodic, which causes them to coherently interefere with each other. If our sampling pattern was aperiod, it would be very unlikely for different signals to have the exact same value at every sample. The easiest way to break coherency is to use a random sampling pattern.
In [7]:
demo.randomized()
The two signals are distinguishable when sampled randomly. This brings us to the second requirement of CS: incoherence.
In [8]:
def random_wrapper(downsample_factor):
demo.random_sample(downsample_factor)
sample_slider = widgets.IntSliderWidget(min=1, max=10, step=1, value=6)
w=widgets.interactive(random_wrapper,downsample_factor=sample_slider)
display(w)
Note the differences in the Fourier spectrum compared to the aliased result we saw earlier. Since we broke the coherency of our sampling pattern, aliased signals now become random noise. We still don't exactly recover our signal, but now that we have removed the ambiguity of aliasing we can try to recover our signal through denoising.
We have seen that random-sampling turns aliased signals into random noise, which we have methods to deal with. Lets look at some possible methods to recover our full signal.
Thresholding
As a simple first step, let's start by thresholding. Here we assume we know something about the amplitude of the signal coefficients.
In [9]:
def random_wrapper(threshold):
demo.random_sample(threshold=threshold, downsample_factor=5)
threshold_slider = widgets.IntSliderWidget(min=0, max=50, step=5, value=20)
w=widgets.interactive(random_wrapper, threshold=threshold_slider)
display(w)
N-largest
In this case we might know before hand how many sinusoids are in our signal. Lets try keeping only the n-largest coefficients.
In [10]:
def random_wrapper(ncoeffs):
demo.random_sample(ncoeffs=ncoeffs, downsample_factor=5)
coeff_slider = widgets.IntSliderWidget(min=1, max=10, step=1, value=2)
w=widgets.interactive(random_wrapper, ncoeffs=coeff_slider)
display(w)
This tutorial demonstrated the basic principles of compressed sensing. We showed that if a signal has a sparse representation and we sample the signal with an incoherent sampling pattern we can fully recover the signal by acquiring far less samples than conventional Nyquist techniques. CS literature has shown drastic sub-sampling under ideal conditions, but engineering requirements in real world applications are still uncertain.
In [43]:
In [ ]: