qgauss.m

非常好的数字处理教程

function [x,y,D,H] = qgauss(N,var,quanttype,disttype)
	sig = sqrt(var);
	sig2 = var;
	
	switch quanttype
		case 'uniform'
			del = 5*sig/N;	% initial quantization region
			e = 10e-4;
			j = 1;

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

						if first
							first = 0;
						end

					end
				end
				y = [-fliplr(q(end-1,:)) q(end-1,:)];
				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,:),var,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)];
				W = 0;
				Dp = distortion(xp,qp,var,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,:),var,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,:),var,disttype);

						if first
							first = 0;
						end

					end
				end
				y = [-fliplr(q(end-1,:)) q(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-10;
			% 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
			switch disttype
				case 'squared error'
					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 and distortion
							W = 0;
							for i=1:N/2

								if x(i+1) == Inf
									y(i) = ( sig*sqrt(2/pi) )* ...
										( exp(-x(i)^2/(2*var)) )/( 1 - erf(x(i)/(sig*sqrt(2))) );
								else
									y(i) = ( sig*sqrt(2/pi) )* ...
										( exp(-x(i)^2/(2*var)) - exp(-x(i+1)^2/(2*var)) )/ ...
										( erf(x(i+1)/(sig*sqrt(2))) - erf(x(i)/(sig*sqrt(2))) );
								end

							end

							% Find new boundaries

							j = j+1;
							D(j+1) = distortion(x,y,var,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;

						% 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

						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) = ( sig*sqrt(2/pi) )* ...
										( exp(-x(i)^2/(2*var)) )/( 1 - erf(x(i)/(sig*sqrt(2))) );
								else
									y(i) = ( sig*sqrt(2/pi) )* ...
										( exp(-x(i)^2/(2*var)) - exp(-x(i+1)^2/(2*var)) )/ ...
										( erf(x(i+1)/(sig*sqrt(2))) - erf(x(i)/(sig*sqrt(2))) );
								end

							end

							% Find new boundaries
							j = j+1;
							D(j+1) = distortion(x,y,var,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'
					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 and distortion
							W = 0;
							for i=1:N/2

								if x(i+1) == Inf
									y(i) = sig*sqrt(2)*erfinv( ...
									  0.5*( 1 + erf(x(i)/(sig*sqrt(2))) ) );
								else
									y(i) = sig*sqrt(2)*erfinv( ...
									  0.5*( erf(x(i+1)/(sig*sqrt(2))) + ...
										erf(x(i)/(sig*sqrt(2))) ) );
								end

							end

							% Find new boundaries

							j = j+1;
							D(j+1) = distortion(x,y,var,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;

						% 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

						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) = sig*sqrt(2)*erfinv( ...
									  0.5*( 1 + erf(x(i)/(sig*sqrt(2))) ) );
								else
									y(i) = sig*sqrt(2)*erfinv( ...
									  0.5*( erf(x(i+1)/(sig*sqrt(2))) + ...
										erf(x(i)/(sig*sqrt(2))) ) );
								end

							end

							% Find new boundaries
							j = j+1;
							D(j+1) = distortion(x,y,var,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;
	sig = sqrt(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
					if x(i+1) ~= Inf
						D = D + sig/(sqrt(2*pi))*( exp(-x(i)^2/(2*var))* ...
							( x(i) - 2*q(i) ) - exp(-x(i+1)^2/(2*var))* ...
							( x(i+1) - 2*q(i) ) ) + 0.5*(var + q(i)^2)*( ...
							erf(x(i+1)/(sig*sqrt(2))) - erf(x(i)/(sig*sqrt(2))) );
					else
						D = D + sig/(sqrt(2*pi))*( exp(-x(i)^2/(2*var))* ...
						( x(i) - 2*q(end) ) ) + 0.5*(var + q(end)^2)* ( ...
						1 - erf(x(i)/(sig*sqrt(2))) );
					end
				end
				if x(i+1) ~= Inf
					% Calculate overload distortion
					D = D + sig/(sqrt(2*pi))*( exp(-x(i+1)^2/(2*var))* ...
						( x(i+1) - 2*q(end) ) ) + 0.5*(var + q(end)^2)* ( ...
						1 - erf(x(i+1)/(sig*sqrt(2))) );
				end
			else
				x = [0 x];
				for i=1:(N+1)/2
					if x(i+1) ~= Inf
						D = D + sig/(sqrt(2*pi))*( exp(-x(i)^2/(2*var))* ...
							( x(i) - 2*q(i) ) - exp(-x(i+1)^2/(2*var))* ...
							( x(i+1) - 2*q(i) ) ) + 0.5*(var + q(i)^2)*( ...
							erf(x(i+1)/(sig*sqrt(2))) - erf(x(i)/(sig*sqrt(2))) );
					else
						D = D + sig/(sqrt(2*pi))*( exp(-x(i)^2/(2*var))* ...
							( x(i) - 2*q(end) ) ) + 0.5*(var + q(end)^2)* ( ...
							1 - erf(x(i)/(sig*sqrt(2))) );
					end
				end
				if x(i+1) ~= Inf
					% Calculate overload distortion
					D = D + sig/(sqrt(2*pi))*( exp(-x(i+1)^2/(2*var))* ...
						( x(i+1) - 2*q(end) ) ) + 0.5*(var + q(end)^2)* ( ...
						1 - erf(x(i+1)/(sig*sqrt(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 + q(i)/2*( 2*erf(q(i)/sqrt(2*var)) - ...
							erf(x(i)/sqrt(2*var)) - erf(x(i+1)/sqrt(2*var)) ) + ...
							1/sqrt(2*pi*var)*( 2*exp(-q(i)^2/(2*var)) - ...
							exp(-x(i)^2/(2*var)) - exp(-x(i+1)^2/(2*var)) );
					else
						D = D + q(i)/2*( 2*erf(q(i)/sqrt(2*var)) - ...
							erf(x(i)/sqrt(2*var)) - 1 ) + ...
							1/sqrt(2*pi*var)*( 2*exp(-q(i)^2/(2*var)) - ...
							exp(-x(i)^2/(2*var)) );
					end
				end
				if x(i+1) ~= Inf
					% Calculate overload distortion
					D = D + q(i)/2*( 2*erf(q(i)/sqrt(2*var)) - ...
							erf(x(i+1)/sqrt(2*var)) - 1 ) + ...
							1/sqrt(2*pi*var)*( 2*exp(-q(i)^2/(2*var)) - ...
							exp(-x(i+1)^2/(2*var)) );
				end
			else
				x = [0 x];
				for i=1:(N+1)/2
					if x(i+1) ~= Inf
						D = D + q(i)/2*( 2*erf(q(i)/sqrt(2*var)) - ...
							erf(x(i)/sqrt(2*var)) - erf(x(i+1)/sqrt(2*var)) ) + ...
							1/sqrt(2*pi*var)*( 2*exp(-q(i)^2/(2*var)) - ...
							exp(-x(i)^2/(2*var)) - exp(-x(i+1)^2/(2*var)) );
					else
						D = D + q(i)/2*( 2*erf(q(i)/sqrt(2*var)) - ...
							erf(x(i)/sqrt(2*var)) - 1 ) + ...
							1/sqrt(2*pi*var)*( 2*exp(-q(i)^2/(2*var)) - ...
							exp(-x(i)^2/(2*var)) );
					end
				end
				if x(i+1) ~= Inf
					% Calculate overload distortion
					D = D + q(i)/2*( 2*erf(q(i)/sqrt(2*var)) - ...
							erf(x(i+1)/sqrt(2*var)) - 1 ) + ...
							1/sqrt(2*pi*var)*( 2*exp(-q(i)^2/(2*var)) - ...
							exp(-x(i+1)^2/(2*var)) );
				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*(erf(x(i+1)/sqrt(2*var)) - erf(x(i)/sqrt(2*var)));
		H = H - p*log2(p);
	end
	H = 2*H;