Implement compressed sensing with FFT

Try to implement compressed sensing with fft. It work fine for dct but give zeros or very small numbers with fft. The code presented here is an example. I have a signal which is sparse in fft... any tips?

%% Generate signal, DCT of signal
n = 4000;    % points in high resolution signal
t = linspace(0, 4000, n); 
x = 0.5*cos(2* 0.008 * pi * t) + 0.5*cos(2* 0.016 * pi * t) + 0.5*cos(2* 0.064 * pi * t) + 0.5*cos(2* 0.256 * pi * t); 
xt = fft(x); % Fourier transformed signal
PSD = abs(real(xt)).^2/((n-1)*1);  % Power spectral density
%% Randomly sample signal
p = 128; % num. random samples, p=n/32
perm = round(rand(p, 1) * n);
y = x(perm); % compressed measurement
%% Solve compressed sensing problem
Psi = dct(eye(n, n));  % build Psi
Theta = Psi(perm, :);  % Measure rows of Psi
%s = cosamp(Theta,y',10,1.e-10,10); % CS via matching pursuit
%% L1-Minimization using CVX
cvx_begin;
variable s(n); 
    minimize( norm(s,1) ); 
    subject to 
        Theta*s == y';
cvx_end;
xrecon = idct(s);      % reconstruct full signal

Thank you!

Answer 1

[Hitesh]

Hi Karoline,

You need to ensure that nature of the signal and the properties of the “FFT” aligns with the compressed signals. The signals you provided is a sum of several cosine waves with different frequencies.

• If these frequencies do not align precisely with the “FFT” bins, the signal won't be sparse in the “FFT” domain. This misalignment causes the energy to spread across multiple “FFT” bins, reducing sparsity.
• For the “FFT” to work well, the frequencies should ideally match the discrete frequencies that the “FFT” can resolve given the signal length.

Compressed sensing relies on the assumption that the signal is sparse in some transform domain. If the signal is not sparse in the “FFT” domain, it is highly recommended to use transform method that better represents the signal's sparsity, such as the “DCT”. That is why it is easier to find the optimal value using the “DCT” rather than the “FFT”. It often leads to an infeasible state like infinity when using “FFT” for compressed sparse signals.

I have revised the code for finding optimal value in compressed signal using “FFT”. You need to ensure following points to avoid reaching the infeasible state:

• Use "SeDumi" solver as it optimized for sparse signal.
• Normalization:  Ensure that the signal normalization step does not lead to division by zero or very small values.
• NaN and Inf Checks: The code checks for NaN or Inf values in the input data to prevent them from causing optimization failures.

Kindly, refer to the revised code:

% Generate signal
n = 4000;    % Number of points in high-resolution signal
t = linspace(0, 4000, n); 
x = 0.5 * cos(2 * 0.008 * pi * t) + 0.5 * cos(2 * 0.016 * pi * t) + ...
    0.5 * cos(2 * 0.064 * pi * t) + 0.5 * cos(2 * 0.256 * pi * t); 

% Fourier Transform of the signal
xt = fft(x); % Fourier transformed signal
PSD = abs(real(xt)).^2 / ((n-1) * 1);  % Power spectral density

% Randomly sample signal
p = 128; % Number of random samples, p = n/32
perm = randperm(n, p);  % Ensure unique and valid indices
y = x(perm); % Compressed measurement

% Normalize measurements
if max(abs(y)) ~= 0
    y = y / max(abs(y));  % Normalize measurements
else
    error('Normalization factor is zero, leading to division by zero.');
end

% Construct measurement matrix using FFT basis
Psi = fft(eye(n));  % FFT basis
Theta = Psi(perm, :);  % Measurement matrix with selected rows

% Check for NaN or Inf values in Theta and y
if any(isnan(Theta(:))) || any(isinf(Theta(:))) || any(isnan(y)) || any(isinf(y))
    error('Inputs contain NaN or Inf values.');
end

% Solve compressed sensing problem using L1 minimization
cvx_begin;
    cvx_precision high;  % Increase precision
    variable s(n); 
    minimize(norm(s, 1)); 
    subject to 
        norm(Theta * s - y') <= 0.01;  % Allow small tolerance
cvx_end;

% Reconstruct the full signal using IFFT
xrecon = ifft(s);

For more information on “FFT”, refer to the MATLAB documentation,

https://www.mathworks.com/help/matlab/ref/fft.html