phaseplane.m

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

function [TFPlane, TFScaleOut] = PhasePlane(c, TFtype, n, qmf, nTFR, TFScale)
% PhasePlane -- partition phase space by rectangular blocks
%  Usage
%    [TFPlane TFScaleOut] = PhasePlane(c, TFtype, n, qmf, nTFR, TFScale)
%  Inputs
%    c		coefficients in the TF dictionary; column vector
%    TFtype     string: type of TF-packets used ('WP','STAT', 'CP')
%    n		signal length
%    qmf	qmf for calculating WP/STAT phase plane location;
%		ignorable for CP
%    nTFR       number of x-points in phase plane image; default = 256
%    TFScale	scale of phase plane image; default maps maximum to 255
%  Outputs
%    TFPlane	phase plane matrix; nTFR by nTFR
%    TFScaleOut scale of phase plane image
%  Side Effects
%    Image of Time-Frequency Plane, with shading showing Heisenberg
%    cells of the atoms making a significant contribution.
%

global COLOR;

if (nargin < 4),          	% Is Wavelet/STAT Packet Phase Plane
	Corrected = 0;    	%    Position-Corrected?
else
	Corrected = 1;
end
if nargin < 5,
	nTFR = 256;
end
if nargin < 6,
	TFScale = 0;
end
titlestr = ' ';

m = length(c);
L = m / n; D = L - 1;
pkt = reshape(c.^ 2, n, L);
TFPlane = zeros(nTFR,nTFR);

if strcmp(TFtype, 'WP'),
	for d = 0:D,
		for b = 0:(2^d-1),
			ylo = b/2^d; yhi = (b+1)/2^d;
			nylo = 1 + floor((nTFR *ylo));
			nyhi = 1 + floor((nTFR *yhi))-1;
			nyhi = max([nyhi nylo]);
			yind = WrapAround(nylo:nyhi,nTFR);
			coeffs = pkt(packet(d,b,n),d+1);
			nbox = n ./ 2^d ; % = length(coeffs);
			if Corrected,
				pos = CalcWPLocation(d,b,0,qmf,n);
			else 
				pos = 0;
			end
			for j=0:(nbox-1),
				jnew = rem((pos/n)*nbox + j, nbox);
				if coeffs(1+j) > 0,
				   nxlo = 1 + floor((nTFR *jnew)/nbox);
				   nxhi = 1 + floor((nTFR *(jnew+1))/nbox)-1;
				   nxhi = max([nxhi nxlo]);
				   xind = WrapAround(nxlo:nxhi,nTFR);
				   TFPlane(yind, xind) = TFPlane(yind,xind) ...
					+ coeffs(1+j);
				end
			end
		end
	end
elseif strcmp(TFtype,'STAT'),
	coarsest = pkt(:,1);
	pkt(1:(L-1), :) = pkt(2:L, :);
	pkt(L) = coarsest;
	for d = 1:L,
		ylo = 1/2^d; yhi = (1+1)/2^d;
		nylo = 1 + floor((nTFR *ylo));
		nyhi = 1 + floor((nTFR *yhi)) - 1;
		nyhi = max([nyhi nylo]);
		yind = WrapAround(nylo:nyhi,nTFR);
		for b = 0:(2^d-1),
			coeffs = pkt(packet(d,b,n),d);
			nbox = n ./ 2^d ; % = length(coeffs);
			if Corrected,
				pos = CalcWPLocation(d,b,0,qmf,n);
			else 
				pos = 0;
			end
			for j=0:(nbox-1),
				jnew = rem((pos/n)*nbox + j, nbox);
				if coeffs(1+j) > 0,
				   nxlo = 1 + floor((nTFR *jnew)/nbox);
				   nxhi = 1 + floor((nTFR *(jnew+1))/nbox) -1;
				   nxhi = max([nxhi nxlo]);
				   xind = WrapAround(nxlo:nxhi,nTFR);
				   TFPlane(yind, xind) = TFPlane(yind,xind) ...
					+ coeffs(1+j);
				end
			end
		end
	end
elseif strcmp(TFtype,'CP'),
	for d = 0:D,
		for b = 0:(2^d-1),
			ylo = b/2^d; yhi = (b+1)/2^d;
			nylo = 1 + floor((nTFR *ylo));
			nyhi = 1 + floor((nTFR *yhi)) - 1;
			nyhi = max([nyhi nylo]);
			yind = WrapAround(nylo:nyhi,nTFR);
			coeffs = pkt(packet(d,b,n),d+1);
			nbox = n ./ 2^d ; % = length(coeffs);
			for j=0:(nbox-1),
				if coeffs(1+j) > 0,
				   nxlo = 1 + floor((nTFR *j)/nbox);
				   nxhi = 1 + floor((nTFR *(j+1))/nbox) - 1;
				   nxhi = max([nxhi nxlo]);
				   xind = WrapAround(nxlo:nxhi,nTFR);
				   TFPlane(xind, yind) = TFPlane(xind,yind) ...
					+ coeffs(1+j);
				end
			end
		   end
	end
end

if TFScale == 0,
	TFMax = max(max(TFPlane));
	TFPlane = 255 * TFPlane / TFMax;
	TFScaleOut = TFMax;
else
	TFPlane = 255 * TFPlane / TFScale;
	TFScaleOut = TFScale;
end

if exist('COLOR'),
	if COLOR,
		colormap(hot(256));
        end
else
	colormap(1-gray(256));
end
image(linspace(0,1,nTFR),linspace(0,1,nTFR),TFPlane); 
axis('xy'); axis([ 0 1 0 1])
titlestr = ['Phase plane: ' titlestr ];
title(titlestr);
xlabel('Time');
ylabel('Frequency');
title(titlestr)