📄 index.html
字号:
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"><html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Compressed Sampling</title> <meta name="generator" content="MATLAB 7.4"> <meta name="date" content="2008-10-19"> <meta name="m-file" content="index"> <LINK REL="stylesheet" HREF="style.css" TYPE="text/css"> </head> <body> <div class="content"> <h1>Compressed Sampling</h1> <introduction> <p>This numerical tour explores compressed/copressive sensing/sampling reconstruction with both greedy matching pursuit and convex L1 regularization. </p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">Installing toolboxes and setting up the path.</a></li> <li><a href="#8">Compressed Sensing Acquisition</a></li> <li><a href="#17">Compressed Sensing Recovery with Orthogonal Matching Pursuit</a></li> <li><a href="#20">Compressed Sensing on Compressible Signals</a></li> </ul> </div> <h2>Installing toolboxes and setting up the path.<a name="1"></a></h2> <p>You need to download the <a href="../toolbox_general.zip">general purpose toolbox</a> and the <a href="../toolbox_signal.zip">signal toolbox</a>. </p> <p>You need to unzip these toolboxes in your working directory, so that you have <tt>toolbox_general/</tt> and <tt>toolbox_signal/</tt> in your directory. </p> <p><b>For Scilab user:</b> you must replace the Matlab comment '%' by its Scilab counterpart '//'. </p> <p><b>Recommandation:</b> You should create a text file named for instance <tt>numericaltour.sce</tt> (in Scilabe) or <tt>numericaltour.m</tt> to write all the Scilab/Matlab command you want to execute. Then, simply run <tt>exec('numericaltour.sce');</tt> (in Scilab) or <tt>numericaltour;</tt> (in Matlab) to run the commands. </p> <p>Execute this line only if you are using Matlab.</p><pre class="codeinput">getd = @(p)path(path,p); <span class="comment">% scilab users must *not* execute this</span></pre><p>Then you can add these toolboxes to the path.</p><pre class="codeinput"><span class="comment">% Add some directories to the path</span>getd(<span class="string">'toolbox_signal/'</span>);getd(<span class="string">'toolbox_general/'</span>);</pre><h2>Compressed Sensing Acquisition<a name="8"></a></h2> <p>Compressed sensing acquisition corresponds to a random projection of a signal <tt>x</tt> on a few linear vectors. For the recovery of <tt>x</tt> to be possible, it is assumed to be sparsely represented in an orthogonal basis. Up to a change of basis, we suppose that the vector <tt>x</tt> itself is sparse. </p> <p>We begin by generating a sparse signal with <tt>s</tt> randomized values. </p><pre class="codeinput"><span class="comment">% dimension</span>n = 512;<span class="comment">% sparsity</span>s = 17;<span class="comment">% location of the diracs</span>set_rand_seeds(123456,123456);sel = randperm(n); sel = sel(1:s);x = zeros(n,1); x(sel)=1;<span class="comment">% % randomization of the signs and values</span>x = x.*sign(randn(n,1)).*(1-.5*rand(n,1));<span class="comment">% save for future use</span>x0 = x;</pre><p>We create a random Gaussian matrix and normalize its columns. This defines the acquisition operator <tt>U</tt>. The columns are thus random point on the unit sphere of <tt>R^p</tt>. </p><pre class="codeinput"><span class="comment">% number of measures</span>p = 100;<span class="comment">% Gaussian matrix</span>U = randn(p,n);<span class="comment">% normalization</span>U = U ./ repmat( sqrt(sum(U.^2)), [p 1] );</pre><p>We perform random measurements without noise.</p><pre class="codeinput">y = U*x;</pre><p>Compressed sensing recovery corresponds to solving the inverse problem <tt>y=U*x</tt>, which is ill posed because <tt>x</tt> is higher dimensional than <tt>x</tt>. The reconstruction can be perform with L1 minimization, which regularizes the problems by using the sparsity of the solution. </p> <p><tt>min_{x1} norm(x1,1) subject to U*x1=y</tt>. </p> <p>Using the homotopy algorithm, we compute the solution path for a large number of different regularization parameter <tt>lambda</tt> values, and retrieve the last step, which corresponds to the minimum L1 norm solution. </p><pre class="codeinput">options.niter = Inf;[X1,lambda_list,sparsity_list] = perform_homotopy(U,y,options);xbp = X1(:,end);</pre><p>We display the solution. You can see that the location of the Dirac is perfectly recovered.</p><pre class="codeinput">clf;subplot(2,1,1);plot_sparse_diracs(x);set_graphic_sizes([], 20);title(<span class="string">'Original Signal'</span>);subplot(2,1,2);plot_sparse_diracs(xbp);set_graphic_sizes([], 20);title(<span class="string">'Recovered by L1 minimization'</span>);</pre><img vspace="5" hspace="5" src="index_01.png"> <p><i>Exercice 1:</i> (the solution is <a href="../private/inverse_compressed_sensing/exo1.m">exo1.m</a>) For sparsity values <tt>s</tt> in [14,34], measure the ratio of vector that are perfectly recovered. To do this, generate many input signal (something like 100) and perform CS recovery on these signals to see wether they are recovered or not. </p><pre class="codeinput">exo1;</pre><img vspace="5" hspace="5" src="index_02.png"> <h2>Compressed Sensing Recovery with Orthogonal Matching Pursuit<a name="17"></a></h2> <p>Instead of using L1 minimization, it is possible to use matching pursuit using the columns of <tt>U</tt> as atoms of a redundant dictionary. </p><pre class="codeinput">x = x0;y = U*x;s = sum(x~=0);</pre><p>We perform <tt>s</tt> steps of orthogonal pursuit. This also recovers perfectly the signal. </p><pre class="codeinput">options.nbr_max_atoms = s;xomp = perform_omp(U,y,options);clf;subplot(2,1,1);plot_sparse_diracs(x);set_graphic_sizes([], 20);title(<span class="string">'Original Signal'</span>);subplot(2,1,2);plot_sparse_diracs(xomp);set_graphic_sizes([], 20);title(<span class="string">'Recovered by OMP'</span>);</pre><img vspace="5" hspace="5" src="index_03.png"> <p><i>Exercice 2:</i> (the solution is <a href="../private/inverse_compressed_sensing/exo2.m">exo2.m</a>) Solve the compressed sensing reconstruction using <tt>s</tt> steps of orthogonal matching pursuits. </p><pre class="codeinput">exo2;</pre><img vspace="5" hspace="5" src="index_04.png"> <img vspace="5" hspace="5" src="index_05.png"> <img vspace="5" hspace="5" src="index_06.png"> <h2>Compressed Sensing on Compressible Signals<a name="20"></a></h2> <p>Exactely <tt>s</tt>-sparse signals are of little practical use because natural signals are compressible rather than exactely sparse. We study here the relative performance of non-linear approximation and compressive sampling approximation on signals with polynomially decaying coefficients. </p> <p>A typical exemple of signal that is compressible in the Dirac basis is a signal with fast decaying coefficients. We flip here the sign, although it is not relevant (because the sensing matrix is random). </p><pre class="codeinput">n = 512;alpha = 1.5;x = (1:n)'.^(-alpha);x(2:2:n) = -x(2:2:n);</pre><p>Random signal are generated from this template by randomizing the coefficients.</p><pre class="codeinput">clf;subplot(2,1,1);plot_sparse_diracs( x(randperm(n)) );subplot(2,1,2);plot_sparse_diracs( x(randperm(n)) );</pre><img vspace="5" hspace="5" src="index_07.png"> <p>The measure of sparsity here is <tt>alpha</tt>, which should be >1 for compressible signal, and in (1/2,1) for not so compressible signal. </p> <p>The non-linear approximation error with <tt>M</tt> coefficients of such a signal decays like <tt>norm(f-f_M)^2=M^{-2*alpha+1}</tt>, so that the slope in log/log plot is <tt>-2*alpha+1</tt>. </p><pre class="codeinput">xs = sort(abs(x));<span class="keyword">if</span> xs(1)>xs(n) xs = reverse(xs);<span class="keyword">end</span>err_nl = reverse(cumsum( xs.^2 ));err_nl = err_nl/err_nl(1);<span class="comment">% display</span>clf;plot( log10(1:n), log10(err_nl) );axis(<span class="string">'tight'</span>);set_label(<span class="string">'log10(M)'</span>, <span class="string">'log10(|f-f_M|^2)'</span>);</pre><img vspace="5" hspace="5" src="index_08.png"> <p><i>Exercice 3:</i> (the solution is <a href="../private/inverse_compressed_sensing/exo3.m">exo3.m</a>) Compute the <tt>alpha</tt> slope for the wavelet coefficients of a natural image (for instance 'lena'). You should perform a linear fit of the slope of ordered coefficients for <tt>M</tt> for instance in the range <tt>[.005 .1]*n^2</tt> (where <tt>n^2</tt> is the number of pixels). Use the function <tt>polyfit</tt> to compute the linear regression. </p><pre class="codeinput">exo3;</pre><img vspace="5" hspace="5" src="index_09.png"> <p>The compressed sensing recovery from a compressible signal from noiseless measurements can be performed by pure L1 minimization. It is however better to use a relaxed L1 minimzation </p> <p><tt>min_{x1} 1/2*norm(y-U*x1)^2 + lambda*norm(x1,1)</tt></p> <p>for a well chosen <tt>lambda</tt>. </p> <p>We compute the full regularization path for all the <tt>lambda</tt>. </p><pre class="codeinput"><span class="comment">% measures</span>y = U*x;<span class="comment">% L1 resolution</span>options.niter = Inf;[X1,lambda_list,sparsity_list] = perform_homotopy(U,y,options);</pre><p>We can display the L2 compressed sensing approximation error as a function of the iteration step. You can see that it is not strictly decaying. A slight regularization somehow decreases the error. </p><pre class="codeinput">m = size(X1,2);err = sum( (X1 - repmat(x, [1 m])).^2 );clf;plot(err / norm(x)); axis([1 m 0 .2]);set_label(<span class="string">'#step'</span>, <span class="string">'Error'</span>);</pre><img vspace="5" hspace="5" src="index_10.png"> <p>The number of samples is <tt>p</tt>. We can compute the number of <tt>q</tt> of sample that generates the same non-linear approximation error. The ratio <tt>p/q>1</tt> gives the compressed sensing sub-optimality constant. The closer to 1, the better. </p><pre class="codeinput"><span class="comment">% the best possible compressed sensing error</span>[e,i] = compute_min(err(:), 1);xcs = X1(:,i);<span class="comment">% compute the equivalent number of NL terms</span>[tmp,q] = compute_min(abs(e-err_nl(:)), 1);disp(strcat([<span class="string">'Compressed sensing sub-optimality factor: p/q='</span> num2str(p/q,2)]));</pre><pre class="codeoutput">Compressed sensing sub-optimality factor: p/q=9.1</pre><p class="footer"><br> Copyright ® 2008 Gabriel Peyre<br></p> </div> <!--##### SOURCE BEGIN #####%% Compressed Sampling% This numerical tour explores compressed/copressive sensing/sampling% reconstruction with both greedy matching pursuit and convex L1% regularization.%% Installing toolboxes and setting up the path.%%% You need to download the % <../toolbox_general.zip general purpose toolbox>% and the <../toolbox_signal.zip signal toolbox>.%%% You need to unzip these toolboxes in your working directory, so% that you have |toolbox_general/| and |toolbox_signal/| in your% directory.%%% *For Scilab user:* you must replace the Matlab comment '%' by its Scilab% counterpart '//'.%%% *Recommandation:* You should create a text file named for instance% |numericaltour.sce| (in Scilabe) or |numericaltour.m| to write all the% Scilab/Matlab command you want to execute. Then, simply run% |exec('numericaltour.sce');| (in Scilab) or |numericaltour;| (in Matlab)% to run the commands. %%% Execute this line only if you are using Matlab.getd = @(p)path(path,p); % scilab users must *not* execute this%%% Then you can add these toolboxes to the path.% Add some directories to the pathgetd('toolbox_signal/');getd('toolbox_general/');⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -