简单程序

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

%% 产生感知矩阵和稀疏表示矩阵

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');

上面都是没有添加噪声的算法结果，我们给信号添加一些噪声，以时域信号为例，

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; % 添加噪声

%% Remaining code not changed