近年来，有关压缩感知的话题越来越多，诸如“基于xx的压缩感知xx”，“压缩感知在xx的应用”，还有类似“基于压缩感知的xxxx”的题目，层出不穷，如果有一条曲线来描述这个爆炸式的流行，那$NumbersOfPapers \propto e^{(t-2004)}$可能会比较合适，其中t是当前的时间。我本人也是从2008年底才开始接触这个“新名词”，认真的研读了一下当时的几篇原创性的文章，例如 E. J. Candès, J. Romberg, and T. Tao, “Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information,” IEEE T. Inform. Theory., vol. 52, no. 2, pp. 489–509, 2006. R. G. Baraniuk, “Compressive Sensing,” in Proc. of Annual Conf. Information Sciences and Systems, CISS’2008, 2007, vol. 24, no. 4, pp. iv–v. E. J. Candès and E. J. Candes, “The Restricted Isometry Property and Its Implications for Compressed Sensing,” Comptes Rendus Mathematique, vol. 346, no. 9–10, pp. 589–592, 2008. 按照我的理解，压缩感知理论，是一个同Nyquist对等的采样框架，只不过两者对信号的假设不同而已（Nyquist假设信号的带宽是有限的，CS假设信号是稀疏的）。从这一点来看，压缩感知其实也就是比Nyquist采样定理对信号的假设更严格而已，因而其结果当然可以用更少的采样数目来恢复原始信号。而在实际自然界中，往往大多数信号确实又满足稀疏性这个条件，例如图像、声音、雷达信号、地震波等等。这就使得，压缩感知可以被广泛的应用的这些领域。 采样过程：$Y=\Phi f$ 恢复过程：$f^* = \Phi^{-1} Y$ 其中，$\Phi: F\to \mathbb{R}^M$, $\Phi^{-1}: \mathbb{R}^M\to F$, $F\subset \mathcal{H}$。 压缩感知 传统Nyquist采样 信号假设 k-稀疏 带限信号 采样约束 矩阵（RIP） 2倍带宽的脉冲采样 恢复算法 L1规则化 低通滤波 如果只是去研究压缩感知的理论框架，的确非常完美，概括下来主要包括三个方面的问题： 信号的稀疏性：强调对信号本身的约束； 感知矩阵的非相干性：强调对采样方法的约束； 恢复（重建）算法。 信号的稀疏性（k）提供了可压缩的可能性，压缩感知的非相干性(RIP)提出了保证了压缩数据可恢复的条件（与稀疏性k、采样数m和信号长度N相关），而恢复算法将压缩数据恢复出来。 这三个方面的问题，相互依赖，而又各自自成体系。 虽然，漂亮的压缩感知理论框架吸引了大多数研究人员的眼球，但是是否在使用这种技术以前，首先关注一下下面的问题： 为什么压缩感知是从MRI图像处理中演变而来？ 除了MRI图像本身采集效率问题以外，是不是跟MRI的成像原理也有关系？如果有关系，那么是否应该在使用压缩感知前，事先考虑一下当前应用使用压缩感知的可行性？难道我们用压缩感知做普通的视频采集也有很大优势？ PS: 个人感觉这压缩感知简直就是为MRI量身定做，我不知道其他领域是否有跟MRI相似的成像原理，或许也可以拿来就用，但是如果用CS做普通的视频、图像、语音等等等等，有点儿像杀鸡用牛刀了，弄巧成拙。
Nuit Blanche 是压缩感知领域影响力最大的一个博客。可惜在国内无法浏览这个博客（除非“翻墙”）。最近博客的博主Igor把2012年所发表的博文汇总成册，以方便下载。详细说明转载如下： The Nuit Blanche Chronicles 2012 are here . The e-book/pdf lists all the entries on Nuit Blanche from Dec 14, 2011 till Dec 21, 2012 .The formatting is crap, it doesn't allow the links to show up (it's a real shame, if somebody has a better way, I am all ears) but it does two useful things: the blog can be read offline on a tablet, a nice thing when you want to unplug and put some perspective on all this The search feature within a pdf is likely much better than Google's.
最近我和本源等合作了一篇论文，已经投稿到了IEEE SPL。这篇文章里，又一个新的BSBL算法推导出来。该算法尽管性能上稍微弱于BSBL-EM,BSBL-BO等，但是速度非常快，所以比较适合大规模的问题。 文章的信息如下： Fast Marginalized Block SBL Algorithm by Benyuan Liu, Zhilin Zhang, Hongqi Fan, Zaiqi Lu, Qiang Fu 下载链接: http://arxiv.org/abs/1211.4909 摘要： The performance of sparse signal recovery can be improved if both sparsity and correlation structure of signals can be exploited. One typical correlation structure is intra-block correlation in block sparse signals. To exploit this structure, a framework, called block sparse Bayesian learning (BSBL) framework, has been proposed recently. Algorithms derived from this framework showed promising performance but their speed is not very fast, which limits their applications. This work derives an efficient algorithm from this framework, using a marginalized likelihood maximization method. Thus it can exploit block sparsity and intra-block correlation of signals. Compared to existing BSBL algorithms, it has close recovery performance to them, but has much faster speed. Therefore, it is more suitable for recovering large scale datasets.
我们的BSBL算法从月份投出去，到现在终于被IEEE Trans. on Signal Processing接收了。论文如下： Zhilin Zhang, Bhaskar. D. Rao, Extension of SBL Algorithms for theRecovery of Block Sparse Signals with Intra-Block Correlation , to appear in IEEE Trans. on Signal Processing 文章可以在arXiv上下载: http://arxiv.org/abs/1201.0862 BSBL以及相关试验代码可以在这里下载: http://dsp.ucsd.edu/~zhilin/BSBL.html 文章摘要如下： We examine the recovery of block sparse signals and extend the framework in two important directions; one by exploiting signals' intra-block correlation and the other by generalizing signals' block structure. We propose two families of algorithms based on the framework of block sparse Bayesian learning (BSBL). One family, directly derived from the BSBL framework, requires to know the block structure. Another family, derived from an expanded BSBL framework, is based on a weaker assumption on the block structure, and can be used in the case when the block structure is completely unknown. Using these algorithms we show that exploiting intra-block correlation is very helpful in improving recovery performance. These algorithms also shed light on how to modify existing algorithms or design new ones to exploit such correlation to improve performance. 下面是相关的应用工作，均发表或者接收在IEEE Trans. on Biomedical Engineering： Zhilin Zhang, Tzyy-Ping Jung, Scott Makeig, Bhaskar D. Rao, Compressed Sensing for Energy-Efficient Wireless Telemonitoring of Non-Invasive Fetal ECG via Block Sparse Bayesian Learning , IEEE Trans. Biomedical Engineering, accepted Zhilin Zhang, Tzyy-Ping Jung, Scott Makeig, Bhaskar D. Rao, Compressed Sensing of EEG for Wireless Telemonitoring with Low Energy Consumption and Inexpensive Hardware , IEEE Trans. Biomedical Engineering, accepted 至此，BSBL的工作就宣告一个段落了。从BSBL算法研究开始到现今，围绕这个框架总共 1，发表/接收了4篇文章到IEEE期刊，分别是IEEE Journal of Selected Topics in Signal Processing, IEEE Trans. on Signal Processing，和IEEE Trans. on Biomedical Engineering 2，发表了若干篇会议文章（CVPR,MICCAI,ICASSP等等) 3，四篇期刊文章在审 4，拿到了工业界5位数的现金奖励（美元）和6位数年薪的工作offer（美元） 5，帮国内的朋友结合某应用领域申请到了30万（RMB）的funding 6，1 US 专利（在审）
http://anony3721.blog.163.com/blog/static/5119742012423105942528/ 压缩感知（Compressed sensing），也被称为压缩采样(Compressive sampling)或稀疏采样(Sparse sampling)，是一种寻找欠定线性系统的稀疏解的技术。如果一个线性方程组未知数的数目超过方程的数目，这个方程组被称为欠定，并且一般而言有无数个解。但是，如果这个欠定系统只有唯一一个稀疏解，那么我们可以利用压缩感知理论和方法来寻找这个解。值得注意的是，不是所有欠定线性方程组都有稀疏解。 压缩感知利用了很多信号中所存在的冗余（换言之，这些信号并非完全是噪声）。具体而言，很多信号都是稀疏的；在适当的表示域中，它们的很多系数都是等于或约等于零。 在信号获取阶段，压缩感知在信号并不稀疏的域对信号进行线性测量。 为了从线性测量中重构出原来的信号，压缩感知求解一个称为L1-范数正则化的最小二乘问题。从理论上可以证明，在某些条件下，这个正则化最小二乘问题的解正是原欠定线性系统的稀疏解 % sparse_in_frequency.m % %This code demonstrate compressive sensing example. In this %example the signal is sparse in frequency domain and random samples %are taken in time domain. close all; clear all; %setup path for the subdirectories of l1magic % path(path, 'C:MATLAB7l1magic-1.1Optimization'); % path(path, 'C:MATLAB7l1magic-1.1Data'); %length of the signal N=1024; %Number of random observations to take K=256; %Discrete frequency of two sinusoids in the input signal k1=29; k2=100; n=0:N-1; %Sparse signal in frequency domain. x=sin(2*pi*(k1/N)*n)+sin(2*pi*(k2/N)*n); % This code demonstrates the compressive sensing using a sparse signal in frequency domain. The signal consists of summation of % two sinusoids of different frequencies in time domain. The signal is sparse in Frequency domain and therefore K random % measurements are taken in time domain. figure; subplot(2,1,1); plot(x) grid on; xlabel('Samples'); ylabel('Amplitude'); title('Original Signal,1024 samples with two different frequency sinsuoids'); xf=fft(x); subplot(2,1,2); plot(abs(xf)) grid on; xlabel('Samples'); ylabel('Amplitude'); title('Frequency domain, 1024 coefficients with 4-non zero coefficients'); %creating dft matrix B=dftmtx(N); Binv=inv(B); % The inverse discrete Fourier transform matrix, Binv,equals CONJ(dftmtx(N))/N. %Taking DFT of the signal xf = B*x.'; %Selecting random rows of the DFT matrix q=randperm(N); %creating measurement matrix A=Binv(q(1:K),:); % 在IDFT矩阵中任选K=256行 %taking random time measurements y=(A*xf); % 对x的fft后的xf(1024-by-1)的数据做IDFT得到256个时域稀疏采样值，通过plot(real(y))和原来的x对比，注意如何在时域中取K=256个采样值 �lculating Initial guess x0=A.'*y; % 注意：待恢复时域信号xprec的DFT值xp的估计初值x0如何给？ y 是时域稀疏采样值 %Running the recovery Algorithm tic xp=l1eq_pd(x0,A, ,y,1e-5); toc figure; subplot(2,1,1) plot(t,x) grid on; xlabel('Time'); ylabel('Amplitude'); title(sprintf('Original Signal, UWB Pulse RF freq=%g GHz',f/1e9)); subplot(2,1,2) plot(t,real(xp),'r') grid on; xlabel('Time'); ylabel('Amplitude'); title(sprintf('Recovered UWB Pulse Signal with %d random samples',K)); %%%%%%%%%%%%%%%%%%飘逸的分割线%%%%%%%%%%%%%%%%%%%%%%%%%% % 用到的函数 % l1eq_pd.m % % Solve % min_x ||x||_1 s.t. Ax = b % % Recast as linear program % min_{x,u} sum(u) s.t. -u = x = u, Ax=b % and use primal-dual interior point method % % Usage: xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter) % % x0 - Nx1 vector, initial point. % % A - Either a handle to a function that takes a N vector and returns a K % vector , or a KxN matrix. If A is a function handle, the algorithm % operates in "largescale" mode, solving the Newton systems via the % Conjugate Gradients algorithm. % % At - Handle to a function that takes a K vector and returns an N vector. % If A is a KxN matrix, At is ignored. % % b - Kx1 vector of observations. % % pdtol - Tolerance for primal-dual algorithm (algorithm terminates if % the duality gap is less than pdtol). % Default = 1e-3. % % pdmaxiter - Maximum number of primal-dual iterations. % Default = 50. % % cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix. % Default = 1e-8. % % cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored % if A is a matrix. % Default = 200. % % Written by: Justin Romberg, Caltech % Email: jrom@acm.caltech.edu % Created: October 2005 % function xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter) largescale = isa(A,'function_handle'); if (nargin 5), pdtol = 1e-3; end if (nargin 6), pdmaxiter = 50; end if (nargin 7), cgtol = 1e-8; end if (nargin 8), cgmaxiter = 200; end N = length(x0); alpha = 0.01; beta = 0.5; mu = 10; gradf0 = ; % starting point --- make sure that it is feasible if (largescale) if (norm(A(x0)-b)/norm(b) cgtol) disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); AAt = @(z) A(At(z)); = cgsolve(AAt, b, cgtol, cgmaxiter, 0); if (cgres 1/2) disp('A*At is ill-conditioned: cannot find starting point'); xp = x0; return; end x0 = At(w); end else if (norm(A*x0-b)/norm(b) cgtol) disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); opts.POSDEF = true; opts.SYM = true; = linsolve(A*A', b, opts); if (hcond 1e-14) disp('A*At is ill-conditioned: cannot find starting point'); xp = x0; return; end x0 = A'*w; end end x = x0; u = (0.95)*abs(x0) + (0.10)*max(abs(x0)); % set up for the first iteration fu1 = x - u; fu2 = -x - u; lamu1 = -1./fu1; lamu2 = -1./fu2; if (largescale) v = -A(lamu1-lamu2); Atv = At(v); rpri = A(x) - b; else v = -A*(lamu1-lamu2); Atv = A'*v; rpri = A*x - b; end sdg = -(fu1'*lamu1 + fu2'*lamu2); tau = mu*2*N/sdg; rcent = - (1/tau); rdual = gradf0 + + ; resnorm = norm( ); pditer = 0; done = (sdg pdtol) | (pditer = pdmaxiter); while (~done) pditer = pditer + 1; w1 = -1/tau*(-1./fu1 + 1./fu2) - Atv; w2 = -1 - 1/tau*(1./fu1 + 1./fu2); w3 = -rpri; sig1 = -lamu1./fu1 - lamu2./fu2; sig2 = lamu1./fu1 - lamu2./fu2; sigx = sig1 - sig2.^2./sig1; if (largescale) w1p = w3 - A(w1./sigx - w2.*sig2./(sigx.*sig1)); h11pfun = @(z) -A(1./sigx.*At(z)); = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0); if (cgres 1/2) disp('Cannot solve system. Returning previous iterate. (See Section 4 of notes for more information.)'); xp = x; return end dx = (w1 - w2.*sig2./sig1 - At(dv))./sigx; Adx = A(dx); Atdv = At(dv); else w1p = -(w3 - A*(w1./sigx - w2.*sig2./(sigx.*sig1))); H11p = A*(sparse(diag(1./sigx))*A'); opts.POSDEF = true; opts.SYM = true; = linsolve(H11p, w1p); if (hcond 1e-14) disp('Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)'); xp = x; return end dx = (w1 - w2.*sig2./sig1 - A'*dv)./sigx; Adx = A*dx; Atdv = A'*dv; end du = (w2 - sig2.*dx)./sig1; dlamu1 = (lamu1./fu1).*(-dx+du) - lamu1 - (1/tau)*1./fu1; dlamu2 = (lamu2./fu2).*(dx+du) - lamu2 - 1/tau*1./fu2; % make sure that the step is feasible: keeps lamu1,lamu2 0, fu1,fu2 0 indp = find(dlamu1 0); indn = find(dlamu2 0); s = min( ); indp = find((dx-du) 0); indn = find((-dx-du) 0); s = (0.99)*min( ); % backtracking line search suffdec = 0; backiter = 0; while (~suffdec) xp = x + s*dx; up = u + s*du; vp = v + s*dv; Atvp = Atv + s*Atdv; lamu1p = lamu1 + s*dlamu1; lamu2p = lamu2 + s*dlamu2; fu1p = xp - up; fu2p = -xp - up; rdp = gradf0 + + ; rcp = - (1/tau); rpp = rpri + s*Adx; suffdec = (norm( ) = (1-alpha*s)*resnorm); s = beta*s; backiter = backiter + 1; if (backiter 32) disp('Stuck backtracking, returning last iterate. (See Section 4 of notes for more information.)') xp = x; return end end % next iteration x = xp; u = up; v = vp; Atv = Atvp; lamu1 = lamu1p; lamu2 = lamu2p; fu1 = fu1p; fu2 = fu2p; % surrogate duality gap sdg = -(fu1'*lamu1 + fu2'*lamu2); tau = mu*2*N/sdg; rpri = rpp; rcent = - (1/tau); rdual = gradf0 + + ; resnorm = norm( ); done = (sdg pdtol) | (pditer = pdmaxiter); disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e, Primal res = %8.3e',... pditer, tau, sum(u), sdg, norm(rdual), norm(rpri))); if (largescale) disp(sprintf(' CG Res = %8.3e, CG Iter = %d', cgres, cgiter)); else disp(sprintf(' H11p condition number = %8.3e', hcond)); end end
Our paper on compressed sensing of EEG for wireless telemonitoring has been accepted by IEEE Trans. on Biomedical Engineering. Here is the summary of this paper: (1) EEG在时域不是稀疏的，在其它变换域（比如离散余弦变换域，小波域）也不是稀疏的（如上图中左下子图所示）； (2) 对于这种非稀疏信号，我们采用BSBL框架中的算法BSBL-BO来恢复； (3) 恢复的质量不仅用常规的性能评价指标来衡量，还采用独立分量分析（ICA）分解来衡量。当恢复的信号有较大误差时，分解出来的独立分量会出现失真。我们试验显示，BSBL恢复的信号的质量能保证独立分量不会出现明显的失真。 (4) 小波压缩与压缩传感的客观比较在压缩传感领域常常被有意或无意地忽视。文中我们用客观冷静的观点评价了两者的优劣。 这篇文章是我们第二篇关于非稀疏信号重建的文章。也是我们脑机接口项目的第一部。不过尽管刚刚被杂志接收，但是里面的算法对我们实验室来说已经过时了。我们现在有更好的算法，并将在下个月投出去。 Anyway, here is the paper: Zhilin Zhang, Tzyy-Ping Jung, Scott Makeig, Bhaskar D. Rao, Compressed Sensing of EEG for Wireless Telemonitoring with Low Energy Consumption and Inexpensive Hardware ,to appear in IEEE Trans. on Biomedical Engineering Abstract: Telemonitoring of electroencephalogram (EEG) through wireless body-area networks is an evolving direction in personalized medicine. Among various constraints in designing such a system, three important constraints are energy consumption, data compression, and device cost. Conventional data compression methodologies, although effective in data compression, consumes significant energy and cannot reduce device cost. Compressed sensing (CS), as an emerging data compression methodology, is promising in catering to these constraints. However, EEG is non-sparse in the time domain and also non-sparse in transformed domains (such as the wavelet domain). Therefore, it is extremely difficult for current CS algorithms to recover EEG with the quality that satisfies the requirements of clinical diagnosis and engineering applications. Recently, Block Sparse Bayesian Learning (BSBL) was proposed as a new method to the CS problem. This study introduces the technique to the telemonitoring of EEG. Experimental results show that its recovery quality is better than state-of-the-art CS algorithms, and sufficient for practical use. These results suggest that BSBL is very promising for telemonitoring of EEG and other non-sparse physiological signals. The code and related data can be downloaded at: https://sites.google.com/site/researchbyzhang/bsbl , or http://dsp.ucsd.edu/~zhilin/BSBL.html
标签: 压缩感知
