bpfig52.m

％ Atomizer Main Directory, Version .802 里面信号含有分解去噪合成过程的代码 ％---------------------------------------

% bpfig52: BP Figure 5.2 -- BP DeNoising Super Resolution
%----------------------------------------------------------------------
%
%          Basis Pursuit is able to achieve super-resolution, while 
% 	   Matching Pursuit sometimes cannot.
%
% Signal:        TwinSine, SNR = 10
%                x = cos(pi * w1 * t) + cos(pi * w2 * t)
% Signal Length: 256
% Dictionary:    4-fold overcomplete discrete cosine dictionary
% Observations:
%    (a) Signal: TwinSine.
%    (b) The inner product b/w the signal and the basis functions.
%        Clearly MP is going to select the basis function with the highest
%        inner product, which is neither w1 nor w2, but in between.
%        Then the residual messes up the whole thing.
%    (c) MP couldn't resolve w1 and w2.
%    (d) The representation from Basis Pursuit has two picks,
%        correponding to w1 and w2.
%
% Use:
%    bpfig52                   uses the current solver.
%    ATOMIZER_ENGINE = 1998;   selects original solver.
%    ATOMIZER_ENGINE = 2001;   selects later    solver.
%    Default is most recent solver.
%----------------------------------------------------------------------

%----------------------------------------------------------------------
%        1998: (S. Chen) Original script for BP paper.
% 09 Apr 2001: (M. Saunders) Choice of solvers implemented.
%----------------------------------------------------------------------

help bpfig52
bpengine;

n     = 256;   fineness = 4;
m     = n * fineness;
t     = ((1:n)' - .5) / n;
const = (2/n) ^ .5;
x1    = const * cos(pi * (125.2  - 1) / fineness * t);
x2    = const * cos(pi * (127.2  - 1) / fineness * t);
%x1   = const * cos(pi * (126    - 1) / fineness * t);
%x2   = const * cos(pi * (128    - 1) / fineness * t);

x     = x1 + x2;
SNR   = 10;
randn('seed', 123456789);
[x y] = NoiseMaker(x, SNR);

%MOFDeNoise
time  = cputime;
[xrecMOF cMOF] = MOFDN(y, 'DCT', fineness, 0, 0);
time1 = cputime - time;

%MPDeNoise
time  = cputime;
[xrecMP cMP] = MPDN(y, 'DCT', fineness, 0, 0);
time3 = cputime - time;

%BPDeNoise
time  = cputime;
switch ATOMIZER_ENGINE
  case 1998
    [xrecBP cBP] = BPDN_Interior (y, 'DCT', fineness, 0, 0, 1e-2, 1e-2, 1e-2);
  case 2001
    [xrecBP cBP] = BPDN_Interior2(y, 'DCT', fineness, 0, 0);
end
time4 = cputime - time;

fprintf('\n')
fprintf('CPU Running Time of MOFDN = %8.4e\n', time1);
fprintf('CPU Running Time of  MPDN = %8.4e\n', time3);
fprintf('CPU Running Time of  BPDN = %8.4e\n', time4);

%--------------------
% Plots
%--------------------
fprintf('\nConstructing figure(%1g) ...\n', FIGURE)
figure(FIGURE);   clf reset;

freqs = (0:(m-1))'/m;

subplot(3,2,1);   plot(t,x);
                  title('(a) TwinSine')

subplot(3,2,3);   plot(t, y);
                  titlestr = sprintf('(b) Noisy TwinSine, SNR = %g', SNR);
                  title(titlestr);

subplot(3,2,5);   plot(freqs, FastDCTAnalysis(x, fineness));
                  title('(c) DCT transform')
                  xlabel('Frequency/Nyquist')

subplot(3,2,2);   plot(freqs, cMOF);
                  title('(d) MOF Coefs')
                  xlabel('Frequency/Nyquist')

subplot(3,2,4);   plot(freqs, cMP);
                  title('(e) MP Coefs')
                  xlabel('Frequency/Nyquist')

subplot(3,2,6);   plot(freqs, cBP);
                  title('(f) BP Coefs')
                  xlabel('Frequency/Nyquist')

subplot(3,2,5);
AXIS = axis;
axis([(101-1)/(fineness*n) (150-1)/(fineness*n) AXIS(3) AXIS(4)]);
AXIS = axis;
X1   = [(125.2-1)/(fineness * n) (125.2-1)/(fineness * n)];
Y1   = [AXIS(3) AXIS(4)];
X2   = [(127.2-1)/(fineness * n) (127.2-1)/(fineness * n)];
Y2   = Y1;
hold on

plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off

subplot(3,2,2);
axis(AXIS);
hold on

plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off

subplot(3,2,4);
axis(AXIS);
hold on

plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off

subplot(3,2,6)
axis(AXIS);
hold on

plot(X1, Y1, ':');
plot(X2, Y2, ':');
hold off