bpfig41.m

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

% bpfig41: BP Figure 4.1 -- Multi-Scale Edges
%----------------------------------------------------------------------
%
%          BP using stationary wavelet dictionaries resembles the 
%          "multi-scale edges" idea of Mallat and Huang.
%
% Signal:        HeaviSine
% Signal Length: 256
% Dictionary:    Stationary Wavelets with D = log2(n) and qmf = (Symmlet, 8)
% Observations:
%    (a) Signal: HeaviSine.
%    (b) BP coeffs displayed in multi-resolution fashion.
%    (c) Stationary Wavelet Transform in multi-resolution display.
%    (d) Mallat's Multi-Scale edges: local maximas of Stationary
%        wavelet transform, in multi-resolution display.
%        Notice that except for the coarsest resolution, BP coefficients
%        correspond to maxima in (d).
%
% Use:
%    bpfig41                   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 bpfig41
bpengine;

n   = 256;
D   = log2(n);
qmf = MakeONFilter('Symmlet', 8);
x   = InputSignal('HeaviSine', n);

%-------------------------------------
% BP
%-------------------------------------
time  = cputime;
switch ATOMIZER_ENGINE
   case 1998
     %c = BP_Interior (x, 'STAT', D-1, qmf, 0, 1e-3, 1e-3, 1e-3);
      c = BP_Interior (x, 'TIWT', D-1, qmf, 0, 1e-3, 1e-3, 1e-3);
   case 2001
      c = BP_Interior2(x, 'TIWT', D-1, qmf, 0 );
end
time4 = cputime - time;

fprintf('\n')
fprintf('CPU Running Time of  BP = %8.4e\n', time4);

pkt1 = TI2Stat(reshape(c, n, D));

%-------------------------------------
% Multi-Scale Edges
%-------------------------------------
pkt2 = FastSTATAnalysis(x, D-1, qmf);
pkt2 = reshape(pkt2', n, D);
pkt3 = pkt2;

for d = 1:(D-1);
   pkt3(:, d+1) = Maxima(abs(pkt2(:, d+1)));
   pkt3(:, d+1) = pkt3(:, d+1) .* sign(pkt2(:, d+1));
end

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

t = (1:n)' / n;
subplot(2,2,1);   plot(t, x);
                  title('(a) Signal: HeaviSine')
subplot(2,2,2);   PlotStatAuto(pkt1);
                  title('(b) Coefs from BP on HeaviSine')
                  xlabel('Position')
                  ylabel('log(resolution)')
subplot(2,2,3);   PlotStatAuto(pkt2);
                  title('(c) Coefs from CWT on HeaviSine')
                  xlabel('Position')
                  ylabel('log(resolution)')
subplot(2,2,4);   PlotStatAuto(pkt3);
                  title('(d) Mutiscale Edges Representation of HeaviSine')
                  xlabel('Position')
                  ylabel('log(resolution)')