qlaplace.m

非常好的数字处理教程

function [x,y,D,H] = qlaplace(N,var,quanttype,disttype)
	sig = sqrt(var);
	sig2 = var;
	l = sqrt(2/var);
	
	switch quanttype
		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,disttype);

				% 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,disttype);

				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,disttype);

						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,disttype);

						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,disttype);

				% 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,disttype);

				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,disttype);

						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,disttype);

						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
			
			switch disttype
				case 'squared error'
					% 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

							end

							% Find new boundaries
							j = j+1;
							D(j+1) = distortion(x,y,sig2,disttype);
							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,disttype);
							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
					
				case 'absolute error'
										% 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) = 1/l*log(2) - ...
										1/l*log(exp(-l*x(i)));
								else
									y(i) = 1/l*log(2) - ...
										1/l*log(exp(-l*x(i+1)) + exp(-l*x(i)));
								end

							end

							% Find new boundaries
							j = j+1;
							D(j+1) = distortion(x,y,sig2,disttype);
							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) = 1/l*log(2) - ...
										1/l*log(exp(-l*x(i)));
								else
									y(i) = 1/l*log(2) - ...
										1/l*log(exp(-l*x(i+1)) + exp(-l*x(i)));
								end

							end

							% Find new boundaries
							j = j+1;
							D(j+1) = distortion(x,y,sig2,disttype);
							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

	end

	H = entropy(x,var);
	
	
function D = distortion(x,q,var,disttype)	
	D = 0;
	l = sqrt(2/var);
	if q(1) == 0
		N = length(q)*2-1;
	else
		N = length(q)*2;
	end
	switch disttype
		case 'squared error'
			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))*( -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))*( -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
		case 'absolute error'
			if rem(N,2) == 0
				% Granular distortion
				for i=1:N/2
					if x(i+1) ~= Inf
						D = D + exp(-l*q(i))/l + 0.5*( ...
							exp(-l*x(i))*(q(i) - x(i) - 1/l) + ...
							exp(-l*x(i+1))*(q(i) - x(i+1) - 1/l) );
					else
						D = D + exp(-l*q(i))/l + ...
							0.5*exp(-l*x(i))*(q(i) - x(i) - 1/l);
					end
				end
				% Calculate Overload Distortion
				if x(i+1) ~= Inf
					D = D + exp(-l*q(i))/l + ...
							0.5*exp(-l*x(i+1))*(q(i) - x(i+1) - 1/l);
				end
				
			else
				x = [0 x];
				for i=1:(N+1)/2
					if x(i+1) ~= Inf
						D = D + exp(-l*q(i))/l + 0.5*( ...
							exp(-l*x(i))*(q(i) - x(i) - 1/l) + ...
							exp(-l*x(i+1))*(q(i) - x(i+1) - 1/l) );
					else
						D = D + exp(-l*q(i))/l + ...
							0.5*exp(-l*x(i))*(q(i) - x(i) - 1/l);
					end
				end
				% Calculate Overload Distortion
				if x(i+1) ~= Inf
					D = D + exp(-l*q(i))/l + ...
							0.5*exp(-l*x(i+1))*(q(i) - x(i+1) - 1/l);
				end
			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;