Implement compressed sensing with FFT

83 次查看(过去 30 天)
Karoline
Karoline 2024-10-23,10:10
编辑: Karoline 2024-10-25,9:07
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!

回答(1 个)

Hitesh
Hitesh 2024-10-24,11:08
编辑:Hitesh 2024-10-24,11:11
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,

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by