clc
clear all
close all
%% https://blog.csdn.net/xiahouzuoxin/article/details/38820925 本程序的网站
% https://blog.csdn.net/HerosOfEarth/article/details/52550213 cvx 的入门
%% 产生原始信号及相关参数
n = 512; % 信号长度
s = 5; % 稀疏度（从下面知道，分时域和频域的情况）
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); % n 个正太随机排列 ，
f = zeros(n,1); % n*1中 找到随机的s 个数值，
f(tmp(1:s)) = randn(s,1); % 产生时域稀疏的信号 %s 个正太分布的数，
end
%% 加噪声 时域稀疏信号不用稀疏表达矩阵， 在频域内进行稀疏表示， 压缩感知进行恢复时，要进行ifft
% 频域稀疏信号则需要先通过稀疏表达矩阵将信号在频域进行稀疏表示，再压缩感知后进行恢复时也要进行反FFT变换到时域
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); % 时域稀疏正交基 , eye 单位矩阵，
end
A = Phi * Psi; % 恢复矩阵 A = Phi * Psi
%上面完成矩阵，
%% 优化方法1：使用CVX工具进行凸优化
addpath('../../cvx-w64/cvx/');
cvx_startup; % 设置cvx环境 cvx_startup 设置cvx环境，
cvx_begin
variable xp(n) complex; % 如果xp是复数，则添加complex，否则不加
minimize (norm(xp, 1));
subject to
A * xp == y;
cvx_end
%% 对比原信号和
if freq_sparse
xp = real(ifft(full(xp))); % 傅里叶逆变换
else
end
plot(f);
hold on
plot(real(xp), 'r.'); % xp 就是要恢复的信号啊，
legend('Original', 'Recovered');
评论0