getnoisepower.m

非常好的数字处理教程

function D = getNoisePower(x,q,var,dist)
D = 0;
switch dist
	case 'uniform'
		V = sqrt(3*var);
		for i=1:length(q)
			D = D + (x(i+1)-q(i))^3 - (x(i)-q(i))^3;
		end
		D = D/(6*V);
			
	case 'gaussian'
		
		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
		
	case 'laplacian'

		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))*( -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
end
D = 2*D;