Tyler Sheffield EE 416 1a. CODE
% EE 416
Exam 2 – Companding (using MATLAB)
Exam 2 "Compander"
- Tyler Sheffield
function y = Compander(x,kind,ce,par); % define the function % x is the array of values mu = 0; % initializations A = 0; C_w = zeros(1,length(x)); % compressor output w = zeros(1,length(x)); % input holder or expander output wp = 1; % set peak if kind == 'm' % will we do mu-law? mu = par; % assign parameter if ce == 'c' % a compressor? w = x; C_w = ( log10(1 + mu .* abs(w/wp) ) .* sign(w) ) ./ log10( 1 + mu ) ; % expression y = C_w; else % default to expander %w = ( wp / mu ) .* sign(x) .* ( exp( abs(x) .* log(1 + mu) / wp) - 1) ; % MATLAB VERSION w = ( wp / mu ) .* sign(x) .* ( ( (1 + mu) .^ abs(x) ) - 1) ; %% HW VERSION - they are the same y = w; end else % default to A-law version A = par; % assign parameter if ce == 'c' % a compressor? w = x; for i=1:1:length(w) % if abs(w(i)/wp) > 1 % C_w(i) = (w(i)/wp) * sign(w(i)); if abs(w(i)/wp) <= (1/A) % compressor expression in first case C_w(i) = ( A * abs(w(i)/wp) * sign(w(i)) ) / ( 1+ log(A) ); else % default case C_w(i) = ( ( ( 1 + log(A * abs(w(i)/wp) ) ) * sign(w(i)) ) ) / ( 1 + log(A) ) ; end end y = C_w; else % default to expander matlab version for i=1:1:length(x) if abs(x(i)) > 1 w(i) = 1*sign(x(i)); elseif abs(x(i)) <= (1/(1+log(A))) % check range w(i) = ( wp * x(i) * ( 1 + log(A))) / A; % expander expression else % default case expander w(i) = exp( ( abs(x(i)) * (1 + log(A))) - 1 ) * wp * sign(x(i)) / A; end end y = w; end end end
1b.
mu-law Compressor (Blue) and A-law Compressor (Red) 1 0.8 0.6 0.4 0.2
output y
0 -0.2 -0.4 -0.6 -0.8 -1 -1.5
-1
-0.5
0 input x
0.5
1
1.5
mu-law Expander (Blue) and A-law Expander (Red) 1 0.8 0.6 0.4 0.2
output y
0 -0.2 -0.4 -0.6 -0.8 -1 -1.5
-1
-0.5
0 input x
0.5
1
1.5
mu-law Compressor (Blue) and A-law Compressor (Red)
0.53
To the right is a zoomed-in image of the above compressor plot showing that the mu-law compressor output gets higher just a tad faster than it’s A-law counterpart. Based on the description of a compressor as a system that boosts the low-level portions of input, the mu-law compressor must be performing greater compression.
0.52
0.51
0.5
output y
0.49
0.48
0.47
0.46
0.02
0.04
0.06 input x
0.08
0.1
0.12
2. CODE
% Exam 2 part 2 - Non-uniform Quantizer function z = Quant(x,kind,par); % define the function % added for correctly finding quantization error N = 2^6; adjustment = 1 - (1/N); % for limiting range of input x x = max(x,-adjustment*ones(1,length(x))); % doing input limiting x = min(x,adjustment*ones(1,length(x))); orig_input = x; from_comp = Compander(x,kind,'c',par); x = from_comp; worry about it plot(orig_input,x,'g'); hold; % pass input through compressor % comes out limited to -1,1 so don't need to
% now run data through quantizer % inputs and initializations N = 2^6; % N-level quantizer, N >= 2 d = 2/N; % this is delta, the step size i = 1; % index into array Q holding quantization output Q = []; % x_forq = 0; adjustment = 1 - (1/N); % for limiting range of x for i=1:1:length(x) x(i) = max(x(i),-adjustment); x(i) = min(x(i),adjustment); % new way of doing input limiting % N even
if mod (N,2) == 0 % step algorithm Q(i) = d * (floor (x(i) / d) + .5); else % N odd Q(i) = d * ( floor ( (x(i) / d) + .5) ); end end Q(i) = Q(i-1); plot(orig_input,Q,'r'); from_exp = Compander(Q,kind,'e',par); z = from_exp; plot(orig_input,z); end
% pass output through expander
PLOTS
mu-law Nonuniform Quantizer Outputs x,y,z 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1
x (G), y (R), z (B)
-0.8
-0.6
-0.4
-0.2
0 input w
0.2
0.4
0.6
0.8
1
A-law Nonuniform Quantizer outputs x,y,z 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1
x (G), y (R), z (B)
-0.8
-0.6
-0.4
-0.2
0 input w
0.2
0.4
0.6
0.8
1
3. See page one of attached engineering paper. 4. See page two of attached engineering paper.
5. CODE
% EE 416 Exam 2 - RndCompander for part 5 function w = RndCompander(kind,par,n,seed); % % % a b v w % define the function
EE 416 Computer Assignment 1 (embedded in Exam 2 code) Random Number Generation initialize boudnaries = 0; = 1; = []; = [];
for i = 1:1:n % generates n samples of the random variable seed = ( mod( (16807 * seed) ,2147483647) ); % perform mod v(i) = a + (b-a) * (seed / 2147483646); % store results in array x if (i == n) % verify visually final value of seed seed; end end if kind == 'm' % mu version mu = par; if mu == 0 % the following statements define the function f_mu f_mu = 2.*(v- .5); else f_mu = (( 1 ./ sqrt ( 1 - ( ((2*mu*(2+mu)) / ((1+mu)^2)) .* abs(v - .5) ) ) ) -1 ) .* sign(v-.5) ./ mu ; end w = f_mu; % assign output else % A version A = par; for i=1:1:length(v) if (abs(v(i)-.5)) <= (1/(3 - (1/(A^2)))) % the following statements define f_A f_A(i) = (1/A) * ( 3 - (1/(A^2)) ) * (v(i) - .5) ; else f_A(i) = ( (1/A) * (sign(v(i) - .5)) ) / ( sqrt(3 - (2 *( 3- (1/(A^2)) ) * abs(v(i) - .5) ) )); end end w = f_A; % assign output end end
6. CODE
% Exam 2 part 2 - Non-uniform Quantizer function z = Quant(x,kind,par); % define the function % added for correctly finding quantization error N = 2^6; adjustment = 1 - (1/N); % for limiting range of input x x = max(x,-adjustment*ones(1,length(x))); % new way of doing input limiting x = min(x,adjustment*ones(1,length(x))); orig_input = x; from_comp = Compander(x,kind,'c',par); x = from_comp; worry about it % plot(orig_input,x,'g'); % hold; % pass input through compressor % comes out limited to -1,1 so don't need to
% now run data through quantizer % inputs and initializations N = 2^6; % N-level quantizer, N >= 2 d = 2/N; % this is delta, the step size i = 1; % index into array Q holding quantization output Q = []; % x_forq = 0; adjustment = 1 - (1/N); % for limiting range of x for i=1:1:length(x) x(i) = max(x(i),-adjustment); x(i) = min(x(i),adjustment); % new way of doing input limiting % N even
if mod (N,2) == 0 % step algorithm Q(i) = d * (floor (x(i) / d) + .5); else % N odd Q(i) = d * ( floor ( (x(i) / d) + .5) ); end end Q(i) = Q(i-1); % plot(orig_input,Q,'r'); from_exp = Compander(Q,kind,'e',par); z = from_exp; %plot(orig_input,z);
% pass output through expander
% Quantization error portion for Exam 2 part 6 eps = z - orig_input; eps2 = eps .* eps; MSE = sum (eps2) / length(eps2) end
CODE (added at end of RndCompander)
% what to run on the end of this (for Exam 2 part 6) x = w; % ee416_ca_04_uniform; % run through uniform quantizer only % y = Quant(w,'m',255); % run for mu-law y = Quant(w,'a',87.6); % run for A-law
CODE (added to end of uniform quantizer)
% Quantization Error portion for Exam 2 part 6 q = y - x; q2 = q .* q; u_MSE = sum(q2) / length(q2)
MSE Table (all values E-6) probability density p sub mu(*): number of samples uniform MSE 10,000 173.4 100,000 173.7 theoretical 81.38 probability density p sub A(*): number of samples uniform MSE 10,000 84.85 100,000 84.17 theoretical 81.38 mu-law MSE 1.502 .5152 .4268 mu-law MSE 2.950 1.404 ----A-law MSE 1.521 .6564 -----A-law MSE 2.934 1.304 1.159
7. Why do the first table’s non-companding uniform quantizer MSEs not agree with the theoretical value while the second table’s do? To answer this question, we must go back to the generation of the random values. We can see the difference that results in the random values generation by examining the CDFs produced from the mu-law and A-law functions f(*):
mu-law (Blue) and A-law (Red) Random Variable CDFs 1 0.9 0.8 0.7 0.6
P(w <= W)
0.5 0.4 0.3 0.2 0.1 0 -1
-0.8
-0.6
-0.4
-0.2
0 w
0.2
0.4
0.6
0.8
1
The mu-law CDF hugs the top-bottom axis closer, meaning that more of the values of w in this case have smaller variation and the samples are tighter. Since N is even and there is really no ‘0’ quantizer input level, this will result in the uniform quantizer having a more difficult time with so many values hovering near 0. This results in a larger MSE.
Session Transcript ( some parts copied from Command History, others from the command line)
x=[-1.2:.0001:1.2]; y = Compander(x,'m','c',255); plot(x,y,'b'); hold y = Compander(x,'a','c',87.6); plot(x,y,'r'); x=[-1.2:.0001:1.2]; y = Quant(x,'m',255); y = Quant(x,'a',87.6); >> clear; >> w = RndCompander('m',255,10000,1); Current plot released MSE = 1.5016e-006 >> w = RndCompander('a',87.6,10000,1); Current plot held MSE = 2.9490e-006 >> w = RndCompander('m',255,100000,1); MSE = 5.1516e-007 >> w = RndCompander('a',87.6,100000,1); MSE = 1.4042e-006 >> w = RndCompander('a',87.6,10000,1); u_MSE = 8.4854e-005 >> w = RndCompander('m',255,10000,1); u_MSE = 1.7338e-004 >> w = RndCompander('m',255,100000,1);
u_MSE = 1.7369e-004 >> w = RndCompander('a',87.6,100000,1); u_MSE = 8.4165e-005 >> w = RndCompander('a',87.6,100000,1); MSE = 1.3044e-006 >> w = RndCompander('a',87.6,10000,1); MSE = 2.9338e-006 >> w = RndCompander('m',255,10000,1); MSE = 1.5208e-006 >> w = RndCompander('m',255,100000,1); MSE = 6.5635e-007