bpfig53.m

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

% bpfig53: BP Figure 5.3 -- TV DeNoising and Dictionary Merger
%----------------------------------------------------------------------
%	
%          BP with Heaviside dictionaries is equivalent to TV DeNoising.
%
% Signal:      Blocks
% Sample Size: 1024
% Dictionary:
%    (1) Heavisides
%    (2) Waves:       Orthogonal Wavelets with qmf=(Symmlet, 8)
%    (3) Waves+Jumps: Orthogonal Wavelets with qmf=(Symmlet, 8)
%                     normalized by TV + Jumps with decay rate .9
% Observations:
%    (a) Signal: Blocks.
%    (b) Noisy Blocks.
%    (c) Plots of sorted coefficients (compression rate in the three
%        different dictionaries).
%    (d) BP Denoising reconstruction using Heavisides (TV DeNoising).
%    (e) BP Denoising reconstruction using Waves.
%    (f) BP Denoising reconstruction using Waves+Jumps.
%
%    Notice that the quality of reconstruction depends on 
%    how well the dictionary compresses the signal.
%    Since Heavisides give the best compression for Blocks,
%    they deliver the best reconstruction.
%
% Use:
%    bpfig53                   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 bpfig53
bpengine;

n    = 1024;   D = log2(n);   par1 = 3;
t    = (1:n)'/n;
qmf1 = MakeONFilter('Symmlet', 8);
qmf2 = MakeONFilter('Haar');
x    = InputSignal('Blocks', n);

%%Noiseless case
%cHS  = BP_Interior(x,'HS',0,0,0,1e-2,1e-2,1e-3);
 cHS  = [x(1); diff(x)];
iz    = cHS == 0;
cHS(iz) = cHS(iz) + 1e-14;
cSymm = IWT_PO(x, par1, qmf1);

dicts = MakeList('POTV', 'JUMP');
par1s = MakeList(par1, 0.9);
par2s = MakeList(qmf1, 0);
par3s = MakeList(0, 0);

switch ATOMIZER_ENGINE
  case 1998
    cMerg = BP_Interior ( x, dicts, par1s, par2s, par3s, 1e-2, 1e-2, 1e-3 );
  case 2001
    cMerg = BP_Interior2( x, dicts, par1s, par2s, par3s );
end

%%DeNoising
randn('seed', 931316785);
SNR   = 7;

%%HeaviSide Dictionary
[x y] = NoiseMaker(x, SNR);

switch ATOMIZER_ENGINE
  case 1998
    [xrecHS crecHS] = BPDN_Interior ( y, 'HS', 0, 0, 0, 1e-1, 1e-1, 1e-3 );
  case 2001
    [xrecHS crecHS] = BPDN_Interior2( y, 'HS', 0, 0, 0 );
end

%%Symmlet
cclean   = SoftThresh( FWT_PO(y, par1, qmf1), (2 * log(n))^.5 ); 
xrecSymm = IWT_PO( cclean, par1, qmf1 );

%%Merge
dicts = MakeList('POTV', 'JUMP');
par1s = MakeList(par1, 0.9);
par2s = MakeList(qmf1, 0);
par3s = MakeList(0, 0);

switch ATOMIZER_ENGINE
  case 1998
    [xrecMerg crecMerg] = ...
       BPDN_Interior ( y, dicts, par1s, par2s, par3s, 1e-1, 1e-1, 1e-3 );
  case 2001
    [xrecMerg crecMerg] = ...
       BPDN_Interior2( y, dicts, par1s, par2s, par3s );
end

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

subplot(3,2,1);   plot(t, x);
                  title('(a) Signal: Blocks');

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

subplot(3,2,4);   plot(t, xrecHS);     % From solve2: BPDN
                  title('(d) BPDeNoise: Heaviside');

subplot(3,2,5);   plot(t, xrecSymm);   % From IWT_PO
                  title('(e) Wavelet Shrinkage: Symmlet');

subplot(3,2,6);   plot(t, xrecMerg);   % From solve3: BPDN
                  title('(f) BPDeNoise: Jump+Wave');

subplot(3,2,1);   AXIS = axis;
subplot(3,2,1);   axis(AXIS);
subplot(3,2,2);   axis(AXIS);
subplot(3,2,4);   axis(AXIS);
subplot(3,2,5);   axis(AXIS);
subplot(3,2,6);   axis(AXIS);
subplot(3,2,3);   semilogy(reverse(sort(abs(cMerg)))/max(abs(cMerg)));
                  axis([0 n 1e-6 1])
                  title('(c) Sorted Coefs');  % From solve1: BP
                  xlabel('Order')
                  ylabel('Amplitude')
                  hold on;

                  semilogy(reverse(sort(abs(cHS  )))/max(abs(cHS  )), ':');
                  semilogy(reverse(sort(abs(cSymm)))/max(abs(cSymm)), '--');
                  MyAxis(100);
                  axes('Position',[.24 .44 .9 .9],'Visible','off');

h1 = text(0, .00, 'Dotted: Heaviside');   set(h1, 'FontSize', 6);
h2 = text(0, .03, 'Dashed: Wave'     );   set(h2, 'FontSize', 6);
h3 = text(0, .06, 'Solid: Jump+Wave' );   set(h3, 'FontSize', 6);
hold off