$25
Here we will attempt to denoise a sparse signal that is corrupted by random noise. Generate a 1×128 vector, x, with 5 non-zero ([1:5]/5) coefficients and permute them randomly. Add random Gaussian noise with standard deviation σ = 0.05 to the signal to obtain y = x + n, where n← σ ∗ randn(1,128). We will solve:
argmin
xˆ
It turns out that this is very easy to solve. Because the variables xˆi, for 1 ≤ i ≤ 128 are independent, we can minimize each of them separately by solving
argmin ,
xˆi
for all 1 ≤ i ≤ 128. Find the closed form solution to each xˆi. Describe what happens when yi is small compared to λ, and when yi is large. Find solutions for λ = {0.01,0.05,0.1,0.2}, and include the plot the error between actual sparse x and the obtained ones.
Random Frequency Domain Sampling and Aliasing
We’ll now explore the importance of incoherent sampling. First, compute the centered Fourier transform of the sparse signal, X = Fx where F is a Fourier transform operator, i.e. X = fftc(x).
In compressed sensing, we undersample the measurements. We measure a subset of k-space, Xu = Fux where Fu is a Fourier transform evaluated only at a subset of frequency domain samples. This is an underdetermined data set for which there is an infinite number of possible signals. However, we do know that the original signal is sparse, so there is hope we will be able to reconstruct it. The theory of compressed sensing suggests random undersampling. To see why, we will look at equispaced undersampling and compare it to random undersampling.
Undersample k-space by taking 32 equispaced samples. Compute the inverse Fourier transform, filling the missing data with zeroes, and multiply by 4 to correct for the fact that we have only 1/4 of the samples, i.e.
Xu = zeros(1,128); Xu(1:4:128) = X(1:4:128); xu = ifftc(Xu)*4 .
Show that this is the Minimum Norm Least Squares solution. Plot the absolute value of the result and comment on it.
Now, undersample k-space by taking 32 samples at random. Compute the zero-filled inverse Fourier transform and multiply by 4 again:
Xr = zeros(1,128); prm = randperm(128); Xr(prm(1:32)) = X(prm(1:32)); xr = ifftc(Xr)*4 .
Plot the absolute value, and comment on it. How does this resemble the denoising problem? Note that the "noise" is not really noise, but incoherent aliasing that is contributed by the signal itself. Therefore, we might be able EXACTLY recover the sparse signal.
Reconstruction from Randomly Sampled Frequency Domain Data
Inspired by the denoising example, we will add an 1-norm penalty and solve
argmin Fuxˆ
xˆ
In this case, xˆ is the estimated signal, Fuxˆ is the undersampled Fourier transform of the estimate, and y are the samples of the Fourier transform that we have acquired. Now the variables are coupled through the Fourier transform, and there is no closed-form solution. Solve this Lasso formulation by using ADMM. Make a plot of error between the true x and Fuxˆ as a function of the iteration number, plotting the result for each of the λ = {0.01,0.05,0.1} . To plot the signal at each iteration use DRAWNOW after the plot command.
1