Compressive Sensing - Sparse in frequency example

Question

I am learning about compressed sensing. I have a question regarding a common MATLAB "sparse in frequency" example that can be find online, for example here and here.

What confuses me in these examples is that they do not take the subset of time samples from the time signal directly, but by inverse transforming the spectrum of the signal. Is this correct? Does this not already involve knowing the frequency spectrum of the signal ? I thought that the goal was to reconstruct the time signal from a subset of time samples, without knowing the exact shape of the spectrum. Or am I missing something?

Answer 1

I think I understand your confusion. In your first example, it would have been more realistic to take

y = x(q)

for example.

A possible explanation for the approach you see here is that authors of such examples would like to ensure they obtain the exact sparse solution they generate the example from to compare against. Therefore they generate the measurements directly from it. It is not "cheating" per se in a mathematical sense, but you can say that it does not reflect the practical reality of taking measurements in the time domain.

What could happen when taking the measurements in the time domain is that if you generate sinusoids with arbitrary frequencies, the signal will appear less sparse in the frequency domain when the frequencies do not match the frequencies of the dictionary. If for example you generate a sinusoid at 1.5 Hz and your DFT has atoms corresponding to the frequencies 1 Hz, 2 Hz, 3 Hz, etc., your signal will not be 1-sparse in this dictionary but will have to use multiple atoms at adjacent frequencies to form the sinusoid at 1.5 Hz. A simple way to avoid this is to generate a signal of the prescribed sparsity directly in the frequency domain and then convert it to the time domain as part of taking the measurements.

Technically, it does not involve already knowing the frequency spectrum of the signal, because X is not known in the reconstruction and the sufficient frequency information in the classic sampling-theoretical sense is lost in the multiplication

y = A*X

because A is a wide matrix.

Answer 2

I think I have solved the problem. The examples are actually not very well written. Below I have modified this version of the example to fit the notation of the paper from Candes and Waking. I am now working to extend it to the case in which we also have a sensing base phi (other then identity, as in the example below).

%% Parameters
% Initialize constants and variables
N = 256; % length of signal (must be even in this example)
M = 64;  % size of the subset
S = 2;   % sparsity in the frequency domain


%% Generate signal with S randomly spread sinusoids
% Note that a real-valued sinusoid has two peaks in the frequency domain
freq = randperm(N/2)-1;
freq = freq(1:S).';
n = 0:N-1;
y = sum(sin(2*pi*freq/N*n).', 2);


%% Fourier representation of the signal
fft1 = fft(eye(N));
x = fft1*y;


%% Perform compressed sensing
% M measurements obtained randomly
subset = randperm(N);
subset = subset(1:M);
y_subset = y(subset); % subset of observed sensed values

% Base for representing the signal. This is the inverse of the Fourier
% transform FFT1, which can be computed as the conjugate of FFT normalized
% by N.
psi = conj(fft1)/N;

% Reconstructed Fourier representation
x_hat = l1eq_pd(randn(N,1), psi(subset,:), [], y_subset, 1e-5);

% Reconstructed signal
y_hat = real(psi*x_hat);