% PCM.m
%
% Neste script é mostrado o efeito do erro de quantização na qualidade do sinal
% com digitalização PCM.
%
% Com base numa amostra de som (.wav, não comprimida) digitalizada,
% preferencialmente a 16 bits/amostra, é feita a quantização com 8, 6, 4 e 2
% bits por amostra. Os sinais quantizados pode ser verificados visualmente, mas
% são também gravados para que possamos identificar o efeito do ruído de 
% quantização na qualidade do som.
%
% Posteriormente o sinal é comprimido segundo a lei A, quantizado a 2 e 4 bits 
% e expandido. Verificamos visulamente o efeito de quantização não uniforme
% causado pela compressão e expansão, e podemos ouvir a melhora na qualidade do
% som.

% André Noll Barreto, 09 ago 2010

close all;

%%%%%%%%%%%%%%
% parâmetros
arquivoSom = 'female1.wav';
us_rate = 8e3;

%%%%%%%%%%%%%%%

[som, fs, bits] = wavread(arquivoSom);

t = (0:length(som)-1) / fs;

% subamostrar
nsamples = length(som);
if mod(nsamples,2)~=0
    nsamples = nsamples - 1;
end

us_factor = fs / us_rate;

som_us = som( 1: us_factor : nsamples);

wavwrite(som_us, us_rate, bits, 'teste_undersample');

%//////////////////////////////////////////////////////////////////////////////
% quantização

mp = max(abs(som));

% quantizar sinal em 8 bits/amostra
nbits = 8;
L = 2^nbits;
som_8bits = round(som / mp * L / 2) * mp / L * 2;
wavwrite(som_8bits, fs, bits, 'teste_8_bits');

% quantizar sinal em 6 bits/amostra
nbits = 6;
L = 2^nbits;
som_6bits =  round(som / mp * L / 2) * mp / L * 2;
wavwrite(som_6bits, fs, bits, 'teste_6_bits');

% quantizar sinal em 4 bits/amostra
nbits = 4;
L = 2^nbits;
som_4bits =  round(som / mp * L / 2) * mp / L * 2;
wavwrite(som_4bits, fs, bits, 'teste_4_bits');

% quantizar sinal em 2 bits/amostra
nbits = 2;
L = 2^nbits;
som_2bits =  round(som / mp * L / 2) * mp / L * 2;
wavwrite(som_2bits, fs, bits, 'teste_2_bits');


figure(1);
plot(t,som)

figure(2);
plot(t,som, t,som_8bits, t, som_6bits, t, som_4bits, t, som_2bits)
legend('16 bits', '8bits', '6 bits', '4 bits', '2 bits');


%%///////////////////////////////////////////////////////////////////////////////
%% compressão
x = .01:.01:1;
mu = 1;
y_mu_0 = 1 / (log(1+mu)) * log (1 + mu*x);

mu = 100;
y_mu_100 = 1 / (log(1+mu)) * log (1 + mu*x);

mu = 1000;
y_mu_1000 = 1 / (log(1+mu)) * log (1 + mu*x);

A = 1;
aux = (x < 1/A) .*  x * A + (x >= 1/A) .* (1 + log(A*x));
y_A_1 = aux / (1 + log(A));

A = 87.6;
aux = (x < 1/A) .*  x * A + (x >= 1/A) .* (1 + log(A*x));
y_A_87 = aux / (1 + log(A));

A = 1000;
aux = (x < 1/A) .*  x * A + (x >= 1/A) .* (1 + log(A*x));
y_A_1000 = aux / (1 + log(A));

figure(3)
plot(x, y_mu_0, x, y_mu_100, x, y_mu_1000, x, y_A_1, x, y_A_87, x, y_A_1000);
legend('\mu=0','\mu=255','\mu=1000', 'A=1', 'A=87,6','A=1000');


% compressão A = 87.6
som_sign = sign(som);
som_amp = abs(som) / mp;
som_amp(som==0) = 1/1000;
A = 87.6;

som_comp = ((som_amp < 1/A) .*  som_amp * A + (som_amp >= 1/A) .* (1 + log(A*som_amp))) / (1 + log(A));
som_comp = som_comp .* som_sign;

wavwrite(som_comp, fs, bits, 'teste_comp');


figure(4)
plot(t, som / mp, t, som_comp);
legend('original','comprimido');

figure(5)
plot(t, som / mp)

% quantizar
nbits = 2;
L = 2^nbits;

som_q =  round(som_comp * L / 2) / L * 2;
som_sign = sign(som_q);
som_q = abs(som_q);
som_2bits_A = ((som_q < 1/(1 + log(A))) .*  som_q * (1 + log(A)) + (som_q >= 1/(1 + log(A))) .* (exp(som_q*(1+log(A))-1))) / A * mp;
som_2bits_A = som_2bits_A .* som_sign;


wavwrite(som_2bits_A, fs, bits, 'teste_aula_2bits');


nbits = 8;
L = 2^nbits;

som_q =  round(som_comp * L / 2) / L * 2;
som_sign = sign(som_q);
som_q = abs(som_q);
som_4bits_A = ((som_q < 1/(1 + log(A))) .*  som_q * (1 + log(A)) + (som_q >= 1/(1 + log(A))) .* (exp(som_q*(1+log(A))-1))) / A * mp;
som_4bits_A = som_4bits_A .* som_sign;

wavwrite(som_4bits_A, fs, bits, 'teste_comp_8bits');

figure(6);
plot(t,som, t,som_4bits_A, t, som_2bits_A)
legend('16 bits', '4bits - lei A', '2 bits - lei A');


figure(7);
subplot(1,2,1)
hist(som, 1000);
title('histograma - voz')
subplot(1,2,2)
hist(som_comp, 1000)
title('histograma - voz com compressão')
