Code covered by the BSD License  

Highlights from
Generate figures for Cleve's Corner on Compressed Sensing

image thumbnail

Generate figures for Cleve's Corner on Compressed Sensing

by

 

21 Jul 2010 (Updated )

Demonstration of compressed sensing.

Cleves_Corner_Compressed_Sensing.m
% Generate the figures for Cleve's Corner.
% "Magic" Reconstruction: Compressed Sensing.
% MathWorks News and Notes, Fall, 2010.
% Use "L1 magic" by Justin Romberg.
% http://www.acm.caltech.edu/l1magic 

% f = signal = tone from "1" key on touch tone phone.
% b = random subsample.

Fs = 40000;
t = (1:Fs/8)'/Fs;
f = (sin(2*pi*697*t) + sin(2*pi*1633*t))/2;
n = length(f);
m = ceil(n/10);
k = randperm(n)';
k = sort(k(1:m));
b = f(k);

% Plot f and b.
% Plot idct(f) = inverse discrete cosine transform.

axf = [0 max(t)/4 -1.2 1.2];
axd = [0 n/8 -10 10];
figure(1);
subplot(2,1,1)
plot(t,f,'b-',t(k),b,'k.')
axis(axf);
set(gca,'xtick',.005:.005:.030,'ytick',-1:1, ...
   'xticklabel',{'.005','.010','.015','.020','.025','.030'})
title('f = signal, b = random sample')
subplot(2,1,2)
plot(idct(f))
axis(axd);
set(gca,'xtick',0:100:600)
title('c = idct(f)')
drawnow

% A = rows of DCT matrix with indices of random sample.

A = zeros(m,n);
for i = 1:m
   ek = zeros(1,n);
   ek(k(i)) = 1;
   A(i,:) = idct(ek);
end

% y = l_2 solution to A*y = b.

y = pinv(A)*b;

% x = l_1 solution to A*x = b.
% Use "L1 magic".

x = l1eq_pd(y,A,A',b,5e-3,32);

% Plot x and dct(x).
% Good comparison with f.

figure(2)
subplot(2,1,1)
plot(x)
axis(axd);
set(gca,'xtick',0:100:600)
title('x = {\it l}_1 solution, A*x = b ')
subplot(2,1,2)
plot(t,dct(x))
axis(axf);
set(gca,'xtick',.005:.005:.030,'ytick',-1:1, ...
   'xticklabel',{'.005','.010','.015','.020','.025','.030'})
title('dct(x)')

% Plot y and dct(y).
% Lousy comparison with f.

figure(3)
subplot(2,1,1)
plot(y)
axis(axd);
set(gca,'xtick',0:100:600)
title('y = {\it l}_2 solution, A*y = b ')
subplot(2,1,2)
plot(t,dct(y))
axis(axf);
set(gca,'xtick',.005:.005:.030,'ytick',-1:1, ...
   'xticklabel',{'.005','.010','.015','.020','.025','.030'})
title('dct(y)')

% Play three sounds.

sound(f,Fs)
pause(1)
sound(dct(x),Fs);
pause(1)
sound(dct(y),Fs)

Contact us