数学描述对于压缩感知模型,其中每个量的维度一定要了解(通过维度的变化来理解压缩感知很有效):
从感知测量中知道:M就是测量的维度(M远小于N)。
压缩感知的原信号恢复问题描述为:
由已知条件:
(1) 测量值y,且,其中e为噪声引入
(2) s为k-Sparse信号(k未知)
求解目标:k尽可能小的稀疏表示信号s及对应的。
用数学形式描述为:
e可以看成告诉随机噪声,e~N(0,δ2)。
即是要求s使s的0范数(非0值的个数)最小,但0范数优化问题是很难求解的,于是一帮大牛就证明求解1范数也能逼近和上面相同的效果,而求解2范数及其更高的范数则结果相差越来越大(有些人在研究介于0范数与1范数之间的范数求解方法)。因此可转化为1范数求解:
由拉格朗日乘子法,上面的最优问题可转化成:
上面的最小值求解问题就可以直接通过凸优化解得结果了。
程序分析先下载CVX或spams工具箱,下面的matlab程序分别使用了时域稀疏信号和频域稀疏信号进行测试,两种不同在于时域稀疏信号不用稀疏表达矩阵(因此稀疏表达矩阵使用单位阵即可),而频域稀疏信号则需要先通过稀疏表达矩阵将信号在频域进行稀疏表示,再压缩感知后进行恢复时也要进行反FFT变换到时域。
关于M的选取:测量数M>=K*log(N/K),K是稀疏度,N信号长度,可以近乎完全重构
clcclear allclose all%% 产生原始信号及相关参数n = 512; % 信号长度s = 25; % 稀疏度(从下面知道,分时域和频域的情况)m = 5*s; % 测量长度 M>=S*log(N/S)freq_sparse = 0; % 信号频域稀疏则为1if 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); % 时域稀疏正交基endA = 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))); % 傅里叶逆变换elseendplot(f+noise);hold onplot(real(xp), 'r.');legend('Original', 'Recovered');
- 设程序中的freq_sparse = 0时,观察时域稀疏信号的恢复结果为
在恢复时域稀疏信号时,所使用的测量矩阵Phi为初始化的随机阵,因为本身时域就稀疏,而算法直接在时域进行恢复,所以稀疏表达矩阵就是单位阵Psi=E。上面的时域稀疏恢复结果显得没有误差是因为没有给原始信号添加噪声。
- 设程序中的freq_sparse = 1时,观察频域稀疏信号的恢复结果为
在恢复时域稀疏信号时,所使用的测量矩阵Phi为初始化的随机阵,因为信号只在频域稀疏,所以稀疏变换矩阵为傅里叶变换基,所以稀疏表达矩阵就是Psi = inv(fft(eye(n,n)))。同理,上面的频域稀疏恢复结果显得没有误差是因为没有给原始信号添加噪声。
- 上面都是没有添加噪声的算法结果,我们给信号添加一些噪声,以时域信号为例,
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); % 产生时域稀疏的信号endnoise = random('norm', 0, 0.01, [n 1]);f = f + noise; % 添加噪声%% Remaining code not changed
从下图的恢复结果看,已经能看到一些恢复误差了,
|