bcsfigure03.m

来自「一种基于压缩感知技术的图像重建程序」· M 代码 · 共 75 行

%-----------------------------------------------------------
% This code generates Figure 3 of the following paper:
% "Bayesian Compressive Sensing" (Preprint, 2007).
% This example is modified from l1qc_example.m, an example
% from l1magic.
% Coded by: Shihao Ji, ECE, Duke University
% last change: Jan. 2, 2007
%-----------------------------------------------------------
% "l1qc_logbarrier.m" and "l1qc_newton.m" from l1-magic are
% required for the BP implementation. They can be found at:
% http://www.acm.caltech.edu/l1magic/
%-----------------------------------------------------------

if exist('l1gc_newton')~=2
    disp(' ')
    disp(' ')
    disp('WARNING: "l1qc_logbarrier.m" and "l1qc_newton.m" from l1-magic')
    disp('webpage are required to run this code. They are available at')
    disp('http://www.acm.caltech.edu/l1magic/.')
    disp(' ')
    disp(' ')
    return
end
clear all
rand('state', 1);
randn('state', 2);
%
N = 512; % signal length
T = 20;  % number of spikes
K = 100; % number of measurements
%
% random +/- 1 signal
x = zeros(N,1);
q = randperm(N);
amp = randn(T,1);
x(q(1:T)) = amp*sqrt(T/sum(amp.^2)); % re-scaled to have the same SNR as in Fig. 2.

% projection matrix
A = randn(K,N);
A = A./repmat(sqrt(sum(A.^2,2)),[1,N]);
% noisy observations
sigma = 0.005;
e = sigma*randn(K,1);
y = A*x + e;

% solve by BP
x0 = A'*inv(A*A')*y;
% take epsilon a little bigger than sigma*sqrt(K)
epsilon =  sigma*sqrt(K)*sqrt(1 + 2*sqrt(2)/sqrt(K));
tic;
x_BP = l1qc_logbarrier(x0, A, [], y, epsilon, 1e-3);
t_BP = toc;
fprintf(1,'BP number of nonzero weights: %d\n',sum(x_BP~=0));

% solve by BCS
initsigma2 = std(y)^2/1e2;
tic;
[weights,used,sigma2,errbars] = BCS_fast_rvm(A,y,initsigma2,1e-8);
t_BCS = toc;
fprintf(1,'BCS number of nonzero weights: %d\n',length(used));
x_BCS = zeros(N,1); err = zeros(N,1);
x_BCS(used) = weights; err(used) = errbars;

% reconstruction error
E_BP = norm(x-x_BP)/norm(x);
E_BCS = norm(x-x_BCS)/norm(x);


subplot(3,1,1); plot(x); axis([1 N -max(abs(x))-0.2 max(abs(x))+0.2]); title(['(a) Original Signal']);
subplot(3,1,2); plot(x_BP); axis([1 N -max(abs(x))-0.2 max(abs(x))+0.2]); title(['(b) Reconstruction with BP, K=' num2str(K)]);
subplot(3,1,3); errorbar(x_BCS,err); axis([1 N -max(abs(x))-0.2 max(abs(x))+0.2]); title(['(c) Reconstruction with BCS, K=' num2str(K)]); box on;

disp(['BP: ||I_hat-I||/||I|| = ' num2str(E_BP) ', time = ' num2str(t_BP) ' secs']);
disp(['BCS: ||I_hat-I||/||I|| = ' num2str(E_BCS) ', time = ' num2str(t_BCS) ' secs']);