bpfig21.m

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

% bpfig21: BP Figure 2.1 -- Time-Frequency Phase Plane
%----------------------------------------------------------------------
%
%                rectangles on phase plane <--> atoms in WP/CP
%
% Signal:        Hydrogen
% Signal Length: 128
% Dictionary:    WP with D = log2(n) and qmf = (symmlet, 8)
% Observations:
%    (a) Hydrogen in frequency domain.
%    (c) Hydrogen in time domain.
%    (b) The rectangle corresponds to the time concentration and
%        frequency concentration of Hydrogen.
%----------------------------------------------------------------------

%----------------------------------------------------------------------
%        1998: (S. Chen) Original script for BP paper.
% 09 Apr 2001: (M. Saunders) Minor editing.
%----------------------------------------------------------------------

help bpfig21
bpengine;

n   = 128;          D = log2(n);
d   = floor(D/2);   b = 2^(d-1) - 1;   k = 2^(D-d-1) - 1;
wp  = MakeWaveletPacket(d,b,k,'Symmlet',8,n);
fwp = abs(fft(wp));
t   = (0:(n-1)) ./n;
wpstr = sprintf('WaveletPacket(%i,%i,%i)',d,b,k);


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

subplot(2,2,2);   qmf = MakeONFilter('Symmlet',8);
%ImageAtomicPhase('WP', [1 d b k ], n, wpstr, 256,qmf);
nTFR = 256;   TFPlane = zeros(nTFR);
f    = k;
nbox = n ./ 2^d;
pos  = CalcWPLocation(d,b,0,qmf,n);
f    = rem((pos/n)*nbox + f, nbox);
ylo  = b/2^d;               yhi = (b +1)/2^d;
xlo  = (f - .5) *(2^d)/n;   xhi = (f - .5 + 1)*(2^d)/n;
nylo = 1 + floor((nTFR *ylo));
nyhi = 1 + floor((nTFR *yhi));
nxlo = 1 + floor((nTFR *xlo));
nxhi = 1 + floor((nTFR *xhi));
xind = WrapAround(nxlo:nxhi, nTFR);
yind = WrapAround(nylo:nyhi, nTFR);

TFPlane(yind, xind) = TFPlane(yind, xind) + 1;
TFMax               = max(max(TFPlane));
TFPlane             = TFPlane .* ( 256./TFMax);

colormap(1 - gray(256))
image(linspace(0,1,nTFR),linspace(0,1,nTFR),TFPlane);
axis('xy'); axis([ 0 1 0 1])
xlabel('Time');
ylabel('Frequency');
title('(b) Phase Plane')
hold on

AXIS = axis;
plot([AXIS(1) AXIS(2)], [ylo ylo], ':')
plot([AXIS(1) AXIS(2)], [yhi yhi], ':')
plot([xlo xlo], [AXIS(3) AXIS(4)], ':')
plot([xhi xhi], [AXIS(3) AXIS(4)], ':')
hold off

subplot(2,2,1);
plot(fwp(1:(n/2)),(1:(n/2))./(n/2))
title('(a) Frequency Domain');
xlabel(['|FFT(' wpstr ')|']);
ylabel('Frequency');
hold on

AXIS = axis;
plot([AXIS(1) AXIS(2)], [ylo ylo], ':')
plot([AXIS(1) AXIS(2)], [yhi yhi], ':')
hold off

subplot(2,2,4);
plot(t,wp)
title('(c) Time Domain');
xlabel('Time');
ylabel(wpstr);
AXIS = axis;
hold on

plot([xlo xlo], [AXIS(3) AXIS(4)], ':')
plot([xhi xhi], [AXIS(3) AXIS(4)], ':')
hold off