qcgauss.m

非常好的数字处理教程

function [x,y,D,H] = qcgauss(N,var,e,type)
	sig = sqrt(var);
	sig2 = var;
	
	switch type
		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);
				
				% 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);

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

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

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

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

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

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

						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
			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);
					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);
					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;
	sig = sqrt(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
			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
	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;