qclaplace.m

非常好的数字处理教程

function [x,y,D,H] = qclaplace(N,var,e,type)
	sig = sqrt(var);
	sig2 = var;
	l = sqrt(2/var);

	switch type
		case 'uniform'
			l = sqrt(2/sig2);
			e = 10e-4;
			j = 1;

			first = 1;		% flag for first time
			del = sqrt(sig2)/N;	% initial quantization region

			% For even # regions
			if rem(N,2) == 0
				x(j,:) = [[0:N/2]*del Inf];
				q(j,:) = x(j,1:end-2)+del/2;

				D(j) = distortion(x(j,:),q(j,:),sig2);

				% Expand a bit to see change in distortion
				xp = x(j,:)*(1+e);
				qp = (xp(2)-xp(1))/2+xp(1:end-2);
				Dp = distortion(xp,qp,sig2);

				first = 1;

				if Dp < D(j)	% Expand regions
					j = j+1;
					D(j) = 0;

					while ( first | D(j) < D(j-1) )
						if ~first
							j = j+1;
						end
						x(j,:) = x(j-1,:)*(1+e);
						q(j,:) = (x(j,2)-x(j,1))/2+x(j,1:end-2);

						D(j) = distortion(x(j,:),q(j,:),sig2);

						if first
							first = 0;
						end

					end
				else	% Contract regions

					j = j+1;
					D(j) = 0;

					while ( first | D(j) < D(j-1) )
						if ~first
							j = j+1;
						end
						x(j,:) = x(j-1,:)*(1-e);
						q(j,:) = (x(j,2)-x(j,1))/2+x(j,1:end-2);

						D(j) = distortion(x(j,:),q(j,:),sig2);

						if first
							first = 0;
						end

					end
				end
				y = [-fliplr(q(end-1,:)) q(end-1,:)];
				%x = [-fliplr(x(end-1,1:end)) x(end-1,2:end)];
				x = [-Inf -fliplr(x(end-1,2:end-2)) x(end-1,1:end-2) Inf];
				D = D(end-1);

			% Odd number of regions
			else 
				x(j,1:(N+1)/2+1) = [[1:(N+1)/2]*del Inf]-del/2;
				q(j,:) = x(j,1:end-1)-del/2;
				D(j) = distortion(x(j,:),q(j,:),sig2);

				% Expand a bit to see change in distortion
				xp = x(j,:)*(1+e);
				qp = [0 (xp(2)-xp(1))/2+xp(1:end-2)];
				Dp = distortion(xp,qp,sig2);

				first = 1;

				if Dp < D(j)	% Expand regions
					j = j+1;
					D(j) = 0;

					while ( first | D(j) < D(j-1) )
						if ~first
							j = j+1;
						end
						x(j,:) = x(j-1,:)*(1+e);
						q(j,:) = [0 (x(j,2)-x(j,1))/2+x(j,1:end-2)];

						D(j) = distortion(x(j,:),q(j,:),sig2);

						if first
							first = 0;
						end

					end
				else	% Contract regions
					j = j+1;
					D(j) = 0;

					while ( first | D(j) < D(j-1) )
						if ~first
							j = j+1;
						end
						x(j,:) = x(j-1,:)*(1-e);
						q(j,:) = [0 (x(j,2)-x(j,1))/2+x(j,1:end-2)];

						D(j) = distortion(x(j,:),q(j,:),sig2);

						if first
							first = 0;
						end

					end
				end
				y = [-fliplr(q(end-1,:)) q(end-1,2:end)];
				%x = [-fliplr(x(end-1,2:end)) x(end-1,2:end)];
				x = [-Inf -fliplr(x(end-1,1:end-2)) x(end-1,1:end-2) Inf];
				D = D(end);
			end

		case 'optimal'
			D = [Inf 10e10];
			j = 1;

			%e = 10e-15;
			% find center of mass for each subset and calculate distortion
			% until difference is less than epsilon
			first = 1;		% flag for first time
			
			% For even # regions
			if rem(N,2) == 0
				x = zeros(1,N/2+1);
				x(1) = 0;
				x(N/2+1) = Inf;

				% initial distribution of ranges, offset by random amount
				for i=2:N/2
					x(i) = 2*sig + (i-1)*(4*sig/N) + (sig/N)*(rand-0.5);
				end

				% Considering N even
				while ( first | abs(D(j)-D(j+1)) > e )

					% Find center of mass
					for i=1:N/2

						if x(i+1) == Inf
							y(i) = ( x(i)*exp(-l*x(i)) )/( exp(-l*x(i)) ) + 1/l;
						else
							y(i) = ( x(i+1)*exp(-l*x(i+1)) - x(i)*exp(-l*x(i)) )/ ...
								( exp(-l*x(i+1)) - exp(-l*x(i)) ) + 1/l;
						end
						x1 = x(i)-y(i);
						x2 = x(i+1)-y(i);

					end
					
					% Find new boundaries
					j = j+1;
					D(j+1) = distortion(x,y,sig2);
					x(N/2+1) = Inf;
					for i=1:N/2-1
						x(i+1) = 0.5*(y(i)+y(i+1));
					end

					first = 0;

				end
				y = [-fliplr(y) y];
				x = [-fliplr(x(1:end)) x(2:end)];
				D = D(end);
			else	% odd number

				x = zeros(1,(N-1)/2+1);
				x((N-1)/2+1) = Inf;

				for i=1:(N-1)/2
					x(i) = sig + (i-1)*(4*sig/N) + (sig/N)*(rand-0.5);
				end

				while ( first | abs(D(j)-D(j+1)) > e )

					% Find center of mass and distortion
					W = 0;
					for i=1:(N-1)/2

						if x(i+1) == Inf
							y(i) = ( x(i)*exp(-l*x(i)) )/( exp(-l*x(i)) ) + 1/l;
						else
							y(i) = ( x(i+1)*exp(-l*x(i+1)) - x(i)*exp(-l*x(i)) )/ ...
								( exp(-l*x(i+1)) - exp(-l*x(i)) ) + 1/l;
						end

					end

					% Find new boundaries
					j = j+1;
					D(j+1) = distortion(x,y,sig2);
					x((N-1)/2+1) = Inf;
					for i=1:(N-1)/2
						if i == 1
							x(i) = 0.5*y(i);
						else
							x(i) = 0.5*(y(i-1)+y(i));
						end
					end

					first = 0;

				end
			y = [-fliplr(y) 0 y];
			x = [-fliplr(x(1:end)) x(1:end)];
			D = D(end);
			end

	end

	H = entropy(x,var);
	
function D = distortion(x,q,var)	
	D = 0;
	l = sqrt(2/var);
	if q(1) == 0
		N = length(q)*2-1;
	else
		N = length(q)*2;
	end
	if rem(N,2) == 0
		% Granular distortion
		for i=1:N/2
			x1 = x(i)-q(i);
			x2 = x(i+1)-q(i);
			if x2 ~= Inf
				D = D + 0.5*( exp(-l*x(i+1))*( -x2^2 - 2*x2/l - 2/l^2 ) - ...
					exp(-l*x(i))*( -x1^2 - 2*x1/l - 2/l^2 ) );
			else
				D = D - 0.5*( exp(-l*x(i+1))*( -x1^2 - 2*x1/l - 2/l^2 ) );
			end
		end
		if x2 ~= Inf
			x1 = x(i+1)-q(i);
			D = D - 0.5*( exp(-l*x(i+1))*( -x1^2 - 2*x1/l - 2/l^2 ) );
		end
		
	else
		x = [0 x];
		for i=1:(N+1)/2
			x1 = x(i)-q(i);
			x2 = x(i+1)-q(i);
			if x2 ~= Inf
				D = D + 0.5*( exp(-l*x(i+1))*( -x2^2 - 2*x2/l - 2/l^2 ) - ...
					exp(-l*x(i))*( -x1^2 - 2*x1/l - 2/l^2 ) );
			else
				D = D - 0.5*( exp(-l*x(i+1))*( -x1^2 - 2*x1/l - 2/l^2 ) );
			end
		end
		if x2 ~= Inf
			x1 = x(i+1)-q(i);
			D = D - 0.5*( exp(-l*x(i+1))*( -x1^2 - 2*x1/l - 2/l^2 ) );
		end
	end
	D = 2*D;
	

	
function H = entropy(x,var)
	H = 0;
	N = length(x)-1;
	if rem(N,2) == 0
		x = x(N/2+1:end);
	else
		x = [0 x((N+1)/2+1:end)];
	end
	
	l = sqrt(2/var);
	
	for i=1:length(x)-1
		p = 0.5*(exp(-l*x(i)) - exp(-l*x(i+1)));
		H = H - p*log2(p);
	end
	H = 2*H;