?? index.html
字號(hào):
<!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>Sparse Spikes Deconvolution with Matching Pursuits</title> <meta name="generator" content="MATLAB 7.4"> <meta name="date" content="2008-10-15"> <meta name="m-file" content="index"> <LINK REL="stylesheet" HREF="style.css" TYPE="text/css"> </head> <body> <div class="content"> <h1>Sparse Spikes Deconvolution with Matching Pursuits</h1> <introduction> <p>This numerical tour explores the use of Matching Pursuits to solve seismic sparse spikes deconvolution.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#3">Installing toolboxes and setting up the path.</a></li> <li><a href="#10">Seismic Wavelets and Seismic Imaging</a></li> <li><a href="#20">Matching Pursuit</a></li> <li><a href="#29">Orthogonal Matching Pursuit</a></li> <li><a href="#34">Exact Recovery Condition</a></li> </ul> </div> <p>Sparse spikes deconvolution is one of the oldest inverse problems, that is a stilized version of recovery in seismic imaging. The ground is modeled as a 1D profile <tt>x</tt>, mostly zeros with a few spikes accounting for interfaces between layers in the ground. The observation is the convolution <tt>D*x</tt> of <tt>x</tt> against a wavelet filter, where <tt>D(:,1)</tt> is the basis wavelet function (<tt>D</tt> is the convolution operator. </p> <p>The goal of sparse spike deconvolution is to recover an approximation of <tt>x</tt> given noisy measurement <tt>y=D*x+w</tt>. Since the convolution destroys many low and high frequencies, this requires some prior information to regularize the inverse problem. Since <tt>x</tt> is composed of a few spikes, the sparsity is a wonderful prior, that is exploited either by L1 minimization (basis pursuit) or by non-linear greedy processes (matching pursuits). </p> <h2>Installing toolboxes and setting up the path.<a name="3"></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>Seismic Wavelets and Seismic Imaging<a name="10"></a></h2> <p>We first see the forard measurement process in action.</p><pre class="codeinput">set_rand_seeds(123456,21);</pre><p>First we load a seismic filter, which is a second derivative of a Gaussian.</p><pre class="codeinput"><span class="comment">% dimension of the signal</span>n = 1024;<span class="comment">% width of the filter</span>s = 13;<span class="comment">% second derivative of Gaussian</span>t = -n/2:n/2-1;h = (1-t.^2/s^2).*exp( -(t.^2)/(2*s^2) );h = h-mean(h);<span class="comment">% normalize it</span>h = h/norm(h);<span class="comment">% recenter the filter for periodic boundary conditions</span>h1 = fftshift(h);</pre><p>We display the filter in space (zoom).</p><pre class="codeinput">k = 100;sel = n/2-k+1:n/2+k;clf;plot(t(sel), h(sel));set_graphic_sizes([], 20);title(<span class="string">'Filter h'</span>);axis(<span class="string">'tight'</span>);</pre><img vspace="5" hspace="5" src="index_01.png"> <p>We display its Fourier transform.</p><pre class="codeinput">hf = fftshift(real(fft(h1)));clf;plot(t(sel), hf(sel));set_graphic_sizes([], 20);title(<span class="string">'FFT(h)'</span>);axis(<span class="string">'tight'</span>);</pre><img vspace="5" hspace="5" src="index_02.png"> <p>The actual number of measurement of the seismic imaging is roughly the number of Fourier frequencies above the noise level <tt>sigma</tt>. </p><pre class="codeinput"><span class="comment">% noise level</span>sigma = .06*max(h);<span class="comment">% how much frequencies are removed</span>q = sum( abs(hf)>sigma );disp(strcat([<span class="string">'Approximate number of measures = '</span> num2str(q) <span class="string">'.'</span>]));</pre><pre class="codeoutput">Approximate number of measures = 106.</pre><p>We compute the filtering matrix. To stabilize the recovery, we sub-sample by a factor of 2 the filtering.</p><pre class="codeinput"><span class="comment">% sub-sampling (distance between wavelets)</span>sub = 2;<span class="comment">% number of atoms in the dictionary</span>p = n/sub;<span class="comment">% the dictionary, with periodic boundary conditions</span>[Y,X] = meshgrid(1:sub:n,1:n);D = reshape( h1(mod(X-Y,n)+1), [n p]);</pre><p>Now we create the ideal signal <tt>x</tt> we would like to recover. We design it so that there is a deacreasing distance between the spikes. It makes the recovery harder and harder from left to right. The amplitude of the spikes and the signs are randoms </p><pre class="codeinput"><span class="comment">% spacing min and max between the spikes.</span>m = 5; M = 40;k = floor( (p+M)*2/(M+m) )-1;spc = linspace(M,m,k)';<span class="comment">% location of the spikes</span>sel = round( cumsum(spc) );sel(sel>p) = [];<span class="comment">% randomization of the signs and values</span>x = zeros(p,1);si = (-1).^(1:length(sel))'; si = si(randperm(length(si)));<span class="comment">% creating of the sparse spikes signal.</span>x(sel) = si;x = x .* (1-rand(p,1)*.5);</pre><p>Display the spikes</p><pre class="codeinput">clf;plot_sparse_diracs(x);title(<span class="string">'Signal x'</span>);</pre><img vspace="5" hspace="5" src="index_03.png"> <p>Sparsity of the signal.</p><pre class="codeinput">disp(strcat([<span class="string">'Sparsity (L0 pseudo-norm) of x = '</span> num2str(sum(x~=0)) <span class="string">'.'</span>]));</pre><pre class="codeoutput">Sparsity (L0 pseudo-norm) of x = 21.</pre><p>Now compute the noisy obervations.</p><pre class="codeinput">w = randn(n,1)*sigma;y = D*x + w;<span class="comment">% display</span>clf;plot(y);set_graphic_sizes([], 20);axis(<span class="string">'tight'</span>);title(<span class="string">'Noisy measurements y=D*x+w'</span>);</pre><img vspace="5" hspace="5" src="index_04.png"> <h2>Matching Pursuit<a name="20"></a></h2> <p>Matching pursuit is a greedy procedure that progressively identify the location of the spikes by looking at atoms that maximaly correlated with the current residual. </p> <p>Initially, the residual is the whole obervations, and the detected signal is zero.</p><pre class="codeinput">R = y;xmp = zeros(p,1);</pre><p>Compute and display the correlation.</p><pre class="codeinput">C = D'*R;clf;plot(abs(C));set_graphic_sizes([], 20);axis(<span class="string">'tight'</span>);title(<span class="string">'|Correlation|'</span>);</pre><img vspace="5" hspace="5" src="index_05.png"> <p>Extract the coefficient with maximal correlation</p><pre class="codeinput">[tmp,I] = compute_max(abs(C));<span class="comment">% update the coefficients</span>xmp(I) = xmp(I) + C(I);<span class="comment">% update the residual</span>R = y-D*xmp;</pre><p>Display the previous and the new residual.</p><pre class="codeinput">clf;subplot(2,1,1);plot(y); axis(<span class="string">'tight'</span>);set_graphic_sizes([], 20);title(<span class="string">'Observation y'</span>);subplot(2,1,2);plot(R); axis(<span class="string">'tight'</span>);set_graphic_sizes([], 20);title(<span class="string">'Residual R'</span>);</pre><img vspace="5" hspace="5" src="index_06.png"> <p><i>Exercice 1:</i> (the solution is <a href="../private/sparsity_seismic_mp/exo1.m">exo1.m</a>) Compute up to M steps of matching pursuit. What do you notice about the locations of the spikes that are well recovered by matching pursuit ? </p><pre class="codeinput">exo1;</pre><img vspace="5" hspace="5" src="index_07.png"> <p>Although several locations of the spikes are well recovered, the values of the spikes are not well estimated. To better recover the values, one needs to find the best values that match the measurements. This requires the least norm solution of the over determined system D(:,sel)*x(sel) = y, where sel is the support recovered by MP. </p><pre class="codeinput">sel = find(xmp~=0);xproj = zeros(p,1);xproj(sel) = D(:,sel) \ y;<span class="comment">% display</span>clf;subplot(2,1,1);plot_sparse_diracs(xmp);title(<span class="string">'Recovered by MP'</span>);subplot(2,1,2);plot_sparse_diracs(xproj);title(<span class="string">'Recovered by backprojected MP'</span>);</pre><img vspace="5" hspace="5" src="index_08.png"> <p><i>Exercice 2:</i> (the solution is <a href="../private/sparsity_seismic_mp/exo2.m">exo2.m</a>) Perform <tt>1.5*M</tt> steps of MP, and at each step compute the back-projection <tt>xproj</tt> of the MP solution <tt>xmp</tt>. Keep the solution <tt>xproj</tt> that minimize the error norm(xproj-x). </p><pre class="codeinput">exo2;</pre><img vspace="5" hspace="5" src="index_09.png"> <p>Display of the solution.</p><pre class="codeinput">err_mp = norm(x-xproj)/norm(x);clf;subplot(2,1,1);plot_sparse_diracs(x);title(<span class="string">'Signal x'</span>);subplot(2,1,2);plot_sparse_diracs(xproj);title([<span class="string">'Recovered by MP, err='</span> num2str(err_mp,3)]);</pre><img vspace="5" hspace="5" src="index_10.png"> <h2>Orthogonal Matching Pursuit<a name="29"></a></h2> <p>Orthogonal matching pursuit improves over matching pursuit by back-projecting at each iteration the matching pursuit solution.</p> <p>The initialization of OMP is the same as with MP.</p><pre class="codeinput">R = y;xomp = zeros(p,1);</pre><p>The coeffcient selection is also done with maximum of correlation.</p><pre class="codeinput">C = D'*R;[tmp,I] = compute_max(abs(C));<span class="comment">% update the coefficients</span>xomp(I) = xomp(I) + C(I);<span class="comment">% perform a back projection of the coefficients</span>sel = find(xomp~=0);xomp = zeros(p,1);xomp(sel) = D(:,sel) \ y;<span class="comment">% update the residual</span>R = y-D*xomp;</pre><p><i>Exercice 3:</i> (the solution is <a href="../private/sparsity_seismic_mp/exo3.m">exo3.m</a>) Implement Orthogonal Matching Pursuit by modifying your implementation of Matching Pursuit. Keep the solution that gives the smallest error <tt>norm(x-xomp)</tt></p><pre class="codeinput">exo3;</pre><img vspace="5" hspace="5" src="index_11.png"> <p>Display of the solution.</p><pre class="codeinput">err_omp = norm(x-xomp)/norm(x);clf;subplot(2,1,1);plot_sparse_diracs(x);title(<span class="string">'Signal x'</span>);subplot(2,1,2);plot_sparse_diracs(xomp);title([<span class="string">'Recovered by OMP, err='</span> num2str(err_omp,3)]);</pre><img vspace="5" hspace="5" src="index_12.png"> <h2>Exact Recovery Condition<a name="34"></a></h2> <p>To study the theoritical properties of pursuits, Joel Tropp introduced a perfect recovery condition (ERC) that measure the conditioning ERC(S) of the support <tt>S=find(x~=0)</tt> of a coefficient. The condition <tt>ERC(S)<1</tt> ensures that any xoefficients supported in <tt>S</tt> is recovered by OMP, up to the noise level (so that the procedure is stable). </p> <p>First we create a signal with well separated Diracs. Such a signal is easy to recover by OMP, and we will check this with the ERC. </p><pre class="codeinput">delta = 30;x1 = zeros(p,1);x1(1:delta:p+1-delta) = 1;</pre><p><tt>ERC(S)</tt> is the maximum L1 norm of the set of inner product between a dual atoms and the set of atoms outside <tt>S</tt>. </p><pre class="codeinput"><span class="comment">% compute the support and the complementary of the support</span>S =find(x1~=0); <span class="comment">% in</span>Sc =find(x1==0); <span class="comment">% out</span><span class="comment">% compute pseudo inverse of atoms within the support</span>D1 = D(:,S);D1 = (D1'*D1)^(-1) * D1';<span class="comment">% compute inner product between dual atoms inside the support</span><span class="comment">% and atoms outside.</span>C = D1 * D(:,Sc);<span class="comment">% compute the maximum L1 norm, which is the ERC</span>ERC = max( sum( abs(C), 1 ) );<span class="comment">% display</span>disp(strcat([<span class="string">'ERC(S)='</span> num2str(ERC,3)]));</pre><pre class="codeoutput">ERC(S)=1.05</pre><p><i>Exercice 4:</i> (the solution is <a href="../private/sparsity_seismic_mp/exo4.m">exo4.m</a>) Compute <tt>ERC(i) = ERC(Si)</tt> for supports <tt>Si</tt> that are separated by an increasing value of <tt>delta</tt>. Check for the minimum <tt>delta</tt> that ensures that <tt>ERC(Si)<1</tt>. </p><pre class="codeinput">exo4;</pre><img vspace="5" hspace="5" src="index_13.png"> <p class="footer"><br> Copyright ® 2008 Gabriel Peyre<br></p> </div> <!--##### SOURCE BEGIN #####%% Sparse Spikes Deconvolution with Matching Pursuits% This numerical tour explores the use of Matching Pursuits to solve% seismic sparse spikes deconvolution.%%% Sparse spikes deconvolution is one of the oldest inverse problems, that% is a stilized version of recovery in seismic imaging.% The ground is modeled as a 1D profile |x|, mostly zeros with a few spikes% accounting for interfaces between layers in the ground. % The observation is the convolution |D*x| of |x| against a wavelet filter,% where |D(:,1)| is the basis wavelet function (|D| is the convolution% operator. %%% The goal of sparse spike deconvolution is to recover an approximation% of |x| given noisy measurement |y=D*x+w|. Since the convolution destroys% many low and high frequencies, this requires some prior information to% regularize the inverse problem. Since |x| is composed of a few spikes,% the sparsity is a wonderful prior, that is exploited either by L1% minimization (basis pursuit) or by non-linear greedy processes (matching% pursuits).%% 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/');%% Seismic Wavelets and Seismic Imaging% We first see the forard measurement process in action.?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -