Compressed Sensing in Julia

Yusri D. Heryanto

In this article, we will explore the application of compressed sensing using Julia. This post is inspired by a blog from Robert Taylor about compressed sensing in Python. Compressed sensing(CS) is a novel sensing paradigm that asserts one can recover signals or images from far fewer samples or measurements than conventional methods use.

CS relies on the principle of sparsity. Most natural signals, such as images and audio, are sparse when the signal is written on an appropriate basis. The commonly used bases for compression in images and audio are the Fourier or wavelet transform.

Mathematically, if a signal inline_formula not implemented is sparse in inline_formula not implemented basis, then instead of measuring y directly (inline_formula not implemented measurement). We can measure a randomly chosen measurement inline_formula not implementedwhere inline_formula not implemented can be dramatically fewer than inline_formula not implemented. The measurements inline_formula not implemented are given by

formula not implemented

The matrix inline_formula not implemented is a sampling matrix. As mentioned above, we assumed inline_formula not implemented is sparse on some appropriate basis inline_formula not implemented:

formula not implemented
blockquote not implemented

Here, the variable inline_formula not implemented is a sparse vector. If the inline_formula not implemented basis is the Fourier basis, inline_formula not implemented can be seen as the inverse Fourier transform of s.

Then, the objective of compressed sensing is to find the sparsest vector inline_formula not implemented from the equation above. This is the optimization problem:

formula not implemented

However, it is non-convex. Fortunately, we can relax the inline_formula not implemented problem to the inline_formula not implemented optimization problem:

formula not implemented

The interested reader can read many papers and articles online such as "A User’s Guide to Compressed Sensing for Communications Systems".

One Dimensional Compressed sensing

First, we explore 1-dimensional compressed sensing. We use the example from Brunton and Kutz. The signal is a mixture of two signals with 2 different frequencies. For the sparse basis choice, we use Discrete Cosine Transform basis.

using Images
using FFTW
using Plots
using LinearAlgebra
using StructuredOptimization
using Distributions
# Create a signal which composed from 2 frequencies
n = 4096 # the number of time points
sampling_rate= n
t = range(0.0, step=1/sampling_rate, length=n)
y = sin.(2* 97 * pi * t) + sin.(2* 777 * pi * t) #freq 194 and 1554
freq = fftfreq(n,sampling_rate);
	p1_ts =plot(t,y,legend=false,title = "Time domain signal")
	p2_ts = plot(t,y,xlim=(0.2,0.3),title = "Inset at [0.2,0.3]")
plot(freq,dct(y),xlim=(0,2000),legend=false,title="Frequency domain signal")

We showed you the signal in the time domain and in the frequency domain. Even though the signal in the time domain is not sparse but the signal in the frequency domain is sparse. This is called the incoherence. The incoherence of the signal is necessary conditions for signal recovery using CS.

# sample a small number of signal
n_sample= div(n,10)
idx = sort(rand(1:n,n_sample)) # sorting is necessary only for nice plotting
#the y_sample here is refer to x in equation above#
y_sample = y[idx];
plot(t[idx],y_sample,seriestype=:scatter,title="Sampling",label="Sampling point",legend=false)

We sample only 10.0 percent of the signal.

Now, using CS theory, we can solve this problem using optimization:

formula not implemented

However, we approached this problem using least absolute shrinkage and selection operator(LASSO) Algorithm:

formula not implemented
s = Variable(n)
~s .= 0.0
lambda = 1e-1
# the function idct(s)[idx] here is same as  performing θ×s
@minimize ls(idct(s)[idx]-y_sample)+lambda*norm(s,1)
((Variable(Float64, (4096,)),), 57)
m_color = palette(:default)[1] #color for recovered signal
t_color = palette(:default)[2] # color for complete signal
plot(plot(freq,dct(y),color=t_color,title="Complete signal (frequency domain)"),
	plot(freq,~s,color=m_color,title="LASSO approach"),

The result is pretty good for the 10.0 percent of data. We can see that the recovery on the boundary is not quite good because violation of the periodic boundary condition requirements of the cosine transform.

Image reconstruction using compressed sensing

Image reconstruction can be seen as 2D or 3D reconstruction. Our example is using an image with three RGB channel colors. We have 2 approaches :(1) First, we can perform the 2D-DCT on each channel or we can perform 3D-DCT. We will show you both approach.

img = load(
)[1:512,1:512] # load image with size 512*512*3
len_img = length(img)
# we will sample about 5% of the image
n_sample_img= round(Int,len_img*0.05);
# We create a mask matrix
mask = zeros(size(img))
mask[sample_idx] .= 1.0
sample_img= mask .* img;
[img Gray.(mask) sample_img]

Here we have the original image, the mask, and the sampled image. We sample about 5.0 percent of the original image.

For the first approach, we extract each channel of the sampled image and flatten it. Then we perform compressed sensing on each channel.

2D Compressed sensing on each channel

#sample the image and flatten it
# Extract each channel of the RGB matrix. This will produced 3*length_of_image matrix. Each row represents each channel.
chanview = convert(Matrix{Float64},channelview(sample_img_array));
r_sample_img_arr= chanview[1,:] # R channel of the sampled image
g_sample_img_arr= chanview[2,:] # G channel of the sampled image
b_sample_img_arr= chanview[3,:]; # B channel of the sampled image
	lambda_2= 3e-2
	r_vx2 = Variable(ny,nx) 
	g_vx2 = Variable(ny,nx) 
	b_vx2 = Variable(ny,nx) 
	@minimize ls(idct(r_vx2)[sample_idx]-r_sample_img_arr)+lambda_2*norm(r_vx2,1) with PANOC()
	@minimize ls(idct(g_vx2)[sample_idx]-g_sample_img_arr)+lambda_2*norm(g_vx2,1) with PANOC()
	@minimize ls(idct(b_vx2)[sample_idx]-b_sample_img_arr)+lambda_2*norm(b_vx2,1) with PANOC()
((Variable(Float64, (512, 512)),), 282)

To see the recovered signal, we inverse transform our variable(i.e. idct(~r_vx2),idct(~g_vx2),idct(~b_vx2))

[colorview(RGB,idct(~r_vx2),idct(~g_vx2),idct(~b_vx2)) img]

3D Compressed sensing on RGB image

For the second approach, we perform 3 Dimensional DCT transform on the 3D array of the image

# Create 3D array from the sampled image
mat = permutedims(channelview(sample_img),(2,3,1));
# Find the linear index of the 3D array where is not zero
sample_idx3=LinearIndices(mat)[findall(x-> x>0, mat)]
# flatten the 3D array
sample_img_arr = mat[sample_idx3];
	lambda_3= 3e-3
	vx3 = Variable(ny,nx,3)
	@minimize ls(idct(vx3)[sample_idx3]-sample_img_arr)+lambda_3*norm(vx3,1) with PANOC()
((Variable(Float64, (512, 512, 3)),), 1000)

The get the image back from vx3, we inverse transform vx3 into n x n x 3 array and then permute the dimensions of array A to show as RGB image.

[colorview(RGB,permutedims(idct(~vx3),(3,1,2))) img]

Here we performed 2 approaches for decompressed sensing of RGB image. The result is recognizable. It shows you the capability of compressed sensing to recover the signal from far less sample.


[1] Brunton, Steven L., and Jose N. Kutz. Data-driven science and engineering: machine learning, dynamical systems, and control. Cambridge, United Kingdom New York, NY: Cambridge University Press, 2019. Print.

[2] Taylor, R. (2016, July 20). Compressed Sensing in Python [Web log post]. Retrieved October 12, 2020, from

[3] Hayashi, K., Nagahara, M., & Tanaka, T. (2013). A User’s Guide to Compressed Sensing for Communications Systems. IEICE Transactions on Communications, E96.B(3), 685–712.

Runtimes (1)