# Compressed Sensing in Julia

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 implementedThe 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 implementedHere, 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 implementedHowever, it is non-convex. Fortunately, we can relax the inline_formula not implemented problem to the inline_formula not implemented optimization problem:

formula not implementedThe 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.

`pkg"add Images FFTW StructuredOptimization Distributions Plots"`

`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);`

`1+1`

`begin`

` 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(p1_ts,p2_ts,legend=false)`

`end`

`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)`

`plot!(t[idx],y_sample,seriestype=:line)`

We sample only 10.0 percent of the signal.

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

formula not implementedHowever, 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`

`ls(idct(s)[idx]-y_sample)+lambda*norm(s,1) `

`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"),`

` xlim=(0,2000),link=:both,legend=false)`

`plot(plot(t,y,color=t_color,title="Complete signal(time domain) "),`

` plot(t,idct(~s,1),color=m_color,title="LASSO approach"),`

` legend=false)`

`plot(plot(t,y,title="Complete signal inset",color=t_color),`

` plot(t,idct(~s,1),title="LASSO approach inset",color=m_color),`

` xlim=(0.2,0.23),legend=false)`

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(bedroom_in_arles.jpg)[1:512,1:512] # load image with size 512*512*3`

`# We create this helper function for random sampling without replacement`

`function rand_wo_replacement(iterable,n_sample)`

` res=Array{eltype(iterable),1}(undef,n_sample)`

` i= 1 `

` while i <= n_sample`

` x = rand(iterable)`

` if x ∈ res`

` continue`

` else `

` res[i] = x`

` end`

` i += 1`

` end`

` return res`

`end;`

`ny,nx=size(img)`

`len_img = length(img)`

`# we will sample about 5% of the image`

`n_sample_img= round(Int,len_img*0.05);`

`sample_idx=rand_wo_replacement(1:len_img,n_sample_img)`

`# 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`

`sample_img_array=img[sample_idx];`

`# 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`

`begin`

` lambda_2= 3e-2`

` r_vx2 = Variable(ny,nx) `

` g_vx2 = Variable(ny,nx) `

` b_vx2 = Variable(ny,nx) `

` `

` ls(idct(r_vx2)[sample_idx]-r_sample_img_arr)+lambda_2*norm(r_vx2,1) with PANOC() `

` `

` ls(idct(g_vx2)[sample_idx]-g_sample_img_arr)+lambda_2*norm(g_vx2,1) with PANOC() `

` `

` ls(idct(b_vx2)[sample_idx]-b_sample_img_arr)+lambda_2*norm(b_vx2,1) with PANOC() `

`end`

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];`

`begin`

` lambda_3= 3e-3`

` vx3 = Variable(ny,nx,3)`

` ls(idct(vx3)[sample_idx3]-sample_img_arr)+lambda_3*norm(vx3,1) with PANOC() `

`end`

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.

## Reference

[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 http://www.pyrunner.com/weblog/2016/05/26/compressed-sensing-python/

[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. https://doi.org/10.1587/transcom.e96.b.685