clc
clear all
close all
%% 产生原始信号及相关参数
n = 512; % 信号长度
s = 25; % 稀疏度（从下面知道，分时域和频域的情况）
m = 5*s; % 测量长度 M>=S*log(N/S)
freq_sparse = 0; % 信号频域稀疏则为1
if freq_sparse
t = [0: n-1]';
f = cos(2*pi/256*t) + sin(2*pi/128*t); % 产生频域稀疏的信号
else
tmp = randperm(n);
f = zeros(n,1);
f(tmp(1:s)) = randn(s,1); % 产生时域稀疏的信号
end
noise = random('norm', 0, 0.01, [n 1]);
f = f + noise; % 添加噪声
%% 产生感知矩阵和稀疏表示矩阵
Phi = sqrt(1/m) * randn(m,n); % 感知矩阵（测量矩阵）
% A = get_A_fourier(n, m);
y = Phi * f; % 通过感知矩阵获得测量值
% Psi 将信号变换到稀疏域
if freq_sparse % 信号频域可以稀疏表示
Psi = inv(fft(eye(n,n))); % 傅里叶正变换，频域稀疏正交基（稀疏表示矩阵）
else % 信号时域可以稀疏表示
Psi = eye(n, n); % 时域稀疏正交基
end
A = Phi * Psi; % 恢复矩阵 A = Phi * Psi
%% 优化方法1：使用CVX工具进行凸优化
% addpath('../../cvx-w64/cvx/');
% cvx_startup; % 设置cvx环境
% cvx_begin
% variable xp(n) complex; % 如果xp是复数，则添加complex，否则不加
% minimize (norm(xp, 1));
% subject to
% A * xp == y;
% cvx_end
%% 优化方法2：使用spams工具箱进行优化
addpath('spams-matlab/build');
param.L = 100;
param.eps = 0.001;
param.numThreads = -1;
A=A./repmat(sqrt(sum(A.^2)),[size(A,1) 1]);
xp = mexOMP(y, A, param); % 正交匹配追踪法Orthogonal Matching Pursuit
%% 对比原信号和
if freq_sparse
xp = real(ifft(full(xp))); % 傅里叶逆变换
else
end
plot(f+noise);
hold on
plot(real(xp), 'r.');
legend('Original', 'Recovered');
