Convolution / Circulant Sampling in 1D and 2D

压缩感知(compressive sensing)在采样前“压缩”信号,从而达到减少采样的数量、时间、成本的目的。理论上,最优的(在order of magnitude和with overwhelming probability意义上)压缩感知方式之一是全随机线性投影。如果用向量x表示一个信号,那么全随机线性投影就是一个矩阵A,它的每个元素是独立的随机数(常见的服从高斯或者伯努利分布)。但是,这样的压缩方式自由度太大,因此很难在机械、物理、电子上简单地实现。然而,压缩感知的目的之一是缩小传感器的尺寸和造价,缩短采样时间。我们又知道,压缩感知必定需要某种全局形式的线性投影。因此,为了在实现中有效的进行压缩采样,我们需要一种物理或者电子上容易实现的线性投影方式,并且要求它依旧有最低的采样要求

Toeplitz/Circulant/卷积形式的随机线性投影比较容易实现。目前,使用它们的具体工作在通信、图像、核磁共振等方面已经出现了。Toeplitz和Circulant矩阵分别是

它们与同样尺寸的高斯随机矩阵相比拥有少得多的自由度,因此也可以说它们更加得不随机。线性投影的随机性在压缩感知很重要,它极大地保证线性投影(encoding linear projection)和信号的稀疏基(sparsifying basis)之间的不一致(incoherent)。这种不一致使得信号可以从很少的采样中被还原。既然Toeplitz和Circulant矩阵比高斯矩阵更加的不随机,那么使用Toeplitz和Circulant矩阵会不会需要更多的采样呢?答案是不会。并且,如果使用了它们,信号还原的算法速度会更快。Toeplitz和Circulant矩阵的产生方法、采样点、最低的采样率、以及快速还原方法在这里有介绍(我的共同作者是Simon Morgan、杨俊峰、张寅)。我们欢迎兴趣的同行索要代码。

最近,一些同行在运行我们的2维代码之后询问我到底什么是2维Circulant Sampling,如何2D图像上实现,跟二维卷积有什么关系,等等。下面,我简单回答一下,并附上Matlab代码作为演示。

This part describes how to implement circulant sampling for both 1D signals and 2D images. It serves an orientation page for algorithms described in this work. First, let us practice generating a standard 1D circulant matrix in Matlab.

Given an n-dimensional 1D vector x and circulant matrix C, one obtains the circulant samples as a matrix-vector product Cx. A circulant matrix is determined by its first row. Other rows are circulant shifts of the first row. In Matlab, one can call C = gallery(‘circul’,v) to generate C whose first row is v. The circulant matrix has a spectral representation: C=(1/n) F^* \mbox{diag}(F^*v) F, or alternatively C=F\mbox{diag}(Fv)F^{-1}, where F and F^* are the Fourier transform and its adjoint, respectively, and v is the first row of C. Therefore, given a vector v, the following Matlab code will yield the same C, C1, and C2:

v = randn(1,10);
n = length(v);
F = dftmtx(n);
C = gallery(‘circul’,v);
C1 = F*diag(F*v(:))/F;
C2 = F’*diag(F’*v(:))*F/n;

To verify, one can run norm(C(:)-C1(:)) and norm(C(:)-C2(:)) and check how small the results are. Often we do not need C but its matrix-vector product Cx. Thanks to the structure of C, having v or its Fourier transform Fv is enough for computing this matrix-vector product. Given x, the following code shall give the same prod1, prod2, prod3, and prod4:

prod1 = C*x;
prod2 = fft(fft(v(:)).*ifft(x));
prod3 = ifft(ifft(v(:)).*fft(x))*n;
temp = conv(v(end:-1:1),x);
prod4 = temp(end-n+1:end);

The formulas for computing prod2 and prod3 are based on the spectral decomposition of C, and they take advantage the fast Fourier transform. prod4 is computed as the convolution between v (in its reversed order) and x. This should not be surprising since Cx consists of the inner products between x and v at all shifted locations. When Cx needs to be computed many times with different x, for better speed I recommend you use the second or third method above and pre-compute fft(v(:)) or ifft(v(:)), respectively.

To move from 1D to 2D, we can extend the 1D spectral representation C=(1/n) F^* \mbox{diag}(F^*v) F by letting F be a 2D Fourier transform and v be a 2D array. Note that C is no longer a matrix but a linear operator on a 2D array, and Cx yield a 2D array consisting of the inner products between x and the 2D array v at its all shifted locations. The shifts are two-way: left-right and up-down. In analogy with the 1D case, Cx becomes a certain 2D convolution for x.

The following code demonstrates five different ways of computing Cx in Matlab, where x and v are randomly generated.

% —— Sec.1 ——-
nRows = ceil(100*rand);
nCols = ceil(100*rand);
x = randn(nRows,nCols);
v = randn(nRows,nCols);
% —— Sec.2 ——-
prod1 = zeros(nRows,nCols);
dot2 = @(x,y) sum(x(:).*y(:));
for jj = 1:nCols
for ii = 1:nRows
prod1(ii,jj) = dot2(v,x);
v = circshift(v,[1 0]);
end
v = circshift(v,[0 1]);
end
% —— Sec.3 ——-
prod2 = fft2(fft2(v).*ifft2(x));
prod3 = ifft2(ifft2(v).*fft2(x))*nRows*nCols;
% —— Sec.4 ——-
psf = flipud(fliplr(circshift(v, [ceil(nRows/2)-1 ceil(nCols/2)-1])));
% —— Sec.5 ——-
prod4 = ifft2(psf2otf(psf).*fft2(x));
prod5 = imfilter(x, psf, ‘conv’, ‘circular’);

Sec. 1 generates x and v of random sizes and endow their entries with random values. They are real valued but they can take complex values too. Cx are computed and stored in prod1, prod2, prod3, prod4, and prod5, all of which have the same size as x and v. They correspond to the moving-mask realization (prod1), FFT realizations (prod2 and prod3), optical realization (prod4), and circular convolution realization (prod5), respectively.

Sec. 2 illustrates how to compute Cx by taking dot-products (i.e., inner products) between x and v at all of its shifted locations. Each inner loop (counted by ii) shifts v downward one step, and each outer loop also shifts v to the right one step. One imagines that x is a static scene and v is a moving mask. Unlike the 1D case, v moves in not just one but both directions. At the end of the nested loops, v returns to its original position as if it is never shifted. Based on the spectral representation of C, Sec. 3 shows how to use fft2 and ifft2 to quickly compute Cx. Sec. 5 provides another two short and fast ways. To use them, however, v must first be converted to a point-spread function (psf). This conversion in Sec. 4 involves shifting v horizontally and vertically about half of its total height and width, respectively, followed by a left-right flip and an upside-down flip. This aligns psf up with Matlab’s psf2otf function and its 2D convolution. The conversion is invertible as v and psf satisfy

v = circshift(fliplr(flipup(psf)), [1-ceil(nRows/2) 1-ceil(nCols/2)]);

Comparing the lines for prod3 and prod4, we also establish ifft2(v)*nRows*nCol = psf2otf(psf). There is nothing magic here really: this explain how Matlab implements psf2otf; it first converts the input psf to v and then applies ifft2(v)*nRows*nCol. If one calls otf = psf2otf(psf,outsize) with outsize larger than the size of psf, then psf must be first padded to that larger size.

Our code RecPC for 2D compressive imaging requires computing Cx, and it uses the form of prod4. The above description hopefully explains how to translate the computation to the physics of random circulant / convolution / Toeplitze compressive sampling and vice versa.

Advertisements
This entry was posted in Compressive Sensing, 稀疏优化. Bookmark the permalink.

5 Responses to Convolution / Circulant Sampling in 1D and 2D

  1. Zhang Cheng says:

    Very Good!

  2. Yousuf says:

    Excellent explanation! Thanks Dr Wotaoyin

  3. YS says:

    Good work, Thanks Dr. WotaoYin

  4. Jiaoyuling says:

    Thank you very much Dr.Yin. And see you in Beijing this summer.

  5. Muhammad says:

    Very informative!!! Thanks

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s