plotphaseplane.sci

小波分解源代码

function PlotPhasePlane(TFType,basis,pkt,titlestr,qmf)
// PlotPhasePlane -- Plot Time-Frequency content of signal
//  Usage
//    PlotPhasePlane(TFType,btree,pkt[,titlestr,qmf]) 
//  Inputs
//    TFType       string 'WP' or 'CP' naming type of analysis
//    btree        basis tree (e.g. as produced by BestBasis)
//    pkt          wavelet or cosine packet table
//    titlestr     optional. signal name for plot title
//    qmf          optional. WP quadrature mirror filter
//
//  Side Effects
//    A plot with rectangles indicating the time-frequency 
//    content of WP or CP components of signal in given basis.
//
//  See Also
//    CPTour, WPTour
//
//  Copyright Aldo I Maalouf  
  
  [lhs,rhs]=argn();
	if (rhs < 5),          // Is Wavelet Packet Phase Plane
		Corrected = 0;        //    Position-Corrected?
	else
		Corrected = 1;
	end

	if rhs < 4,
	   titlestr = ' ';
	end

	[n,L] = size(pkt);
//
// 	Initialize tree traversal stack
//
	dstack = zeros(1,2^L);
	bstack = zeros(1,2^L);
	k = 1;
	dstack(k) = 0;
	bstack(k) = 0;
//
	ss = norm(pkt(:,1));
//
//	Initialize plotting window
//
	LockAxes([ 0 1 0 1]);
	titlestr = ['Phase plane: ' titlestr ];
	xtitle(titlestr);
//
	while(k > 0),
		d = dstack(k); b = bstack(k); k=k-1;
		if(basis(node(d,b)) == 0) ,  // terminal node
			ylo = b/2^d; yhi = (b+1)/2^d;
			coeffs = n .* (pkt(packet(d,b,n),d+1) ./ss) .^2;
			nbox = n ./ 2^d ; 
			if string(TFType)=='WP',
				if Corrected,
					pos = CalcWPLocation(d,b,0,qmf,n);
				else 
					pos = 0;
				end
				for j=0:(nbox-1),
					jpos = rem((pos/n)*nbox + j, nbox);
					DrawHeisenberg(jpos/nbox, (jpos+1)/nbox, ylo,yhi, coeffs(1+j))
				end
			else
				for j=0:(nbox-1),
					DrawHeisenberg( ylo,yhi, j/nbox, (j+1)/nbox, coeffs(1+j))
				end
			end
		else
		   k = k+1;
		   dstack(k) = d+1; bstack(k) = 2*b;
		   k = k+1;
		   dstack(k) = d+1; bstack(k) = 2*b+1;
		 end
	end
	xtitle('','Time','Frequency')
	
	xset("auto clear","on");//UnlockAxes;
	endfunction