Matlab codes

1D test

a test of 1D signal recovery with synthetic dictionaries and a randomly generated sparse signal

2D test

a test that illustrates how to apply a 2D circulant operator to image patches and recover a 2D image from circulat measurements

Fullsize recovery

a test that demonstrates image recovery using a fullsize sensing operator and dictionary; see fullsize recovery

Download source code

How to simulate 1D/2D circulant samples in Matlab

We describe how to simulate circulant sampling on one or two dimensional signals in Matlab. 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 -dimensional 1D vector and circulant matrix , one obtains the circulant samples as a matrix-vector product . 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 whose first row is . The circulant matrix has a spectral representation: , or alternatively , where and are the Fourier transform and its adjoint, respectively, and is the first row of . Therefore, given a vector v, the following Matlab code will produce 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 but its matrix-vector product . Thanks to the structure of , having or its Fourier transform is enough for computing this matrix-vector product. Given , 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 , and they take advantage the fast Fourier transform. prod4 is computed as the convolution between (in its reversed order) and . This should not be surprising since consists of the inner products between and at all shifted locations. When needs to be computed many times with different , for better speed I recommend you use the second or third method above with pre-computed fft(v(:)) or ifft(v(:)), respectively.

Moving from 1D to 2D, we can extend the 1D spectral representation by letting be a 2D Fourier transform and be a 2D array. Note that is no longer a matrix but a linear operator on a 2D array, and yield a 2D array consisting of the inner products between and the 2D array at its all shifted locations. The shifts are two-way: left-right and up-down. In analogy with the 1D case, becomes a certain 2D convolution for .

The following code demonstrates five different ways of computing in Matlab, where and 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 and of random sizes and endow their entries with random values. They are real valued but they can take complex values too. are computed and stored in prod1, prod2, prod3, prod4, and prod5, all of which have the same size as and . 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 by taking dot-products (i.e., inner products) between and at all of its shifted locations. Each inner loop (counted by ii) shifts downward one step, and each outer loop also shifts to the right one step. One imagines that is a static scene and is a moving mask. Unlike the 1D case, moves in not just one but both directions. At the end of the nested loops, returns to its original position as if it is never shifted. Based on the spectral representation of , Sec. 3 shows how to use fft2 and ifft2 to quickly compute . Sec. 5 provides another two short and fast ways. To use them, however, must first be converted to a point-spread function (psf). This conversion in Sec. 4 involves shifting 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 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)nRowsnCol = psf2otf(psf). There is nothing magic here really: this explain how Matlab implements psf2otf; it first converts the input psf to and then applies ifft2(v)nRowsnCol. 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.