% este script simula a amostragem de um sinal de voz e sua reconstrução
%
% Universidade de Brasília
% André Noll Barreto, 11 set 2014


close all

%%%%%%%
% parâmetros
arquivoSom1 = 'female1.wav';
fmin = 16000; % frequência máxima observável na DFT
downsampling_rate = 16;


%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
% super amostrar sinal para aumentar faixa de frequência observável

os = ceil(2*fmin/fs); % fator de superamostragem


t1 = (0:length(som)-1) / fs;
t1_os = (0 : 1/os : (length(som)-1)) / fs; % instantes de amostra
t1_s = 1 / fs / os;


som_os = interp1(t1, som, t1_os);  % som super-amostrado

% estimar densidade espectral de potência
[ som_psd, f] = pwelch( som_os, [], [],[], fs*os);

figure(1);
plot( f, abs(som_psd) );
title('Densidade Espectral de Potência - Mensagem de Voz');
xlabel('frequência (Hz)');

figure(2);
maxdB1 = 10*log10(max(abs(som_psd)));
plot(f, 10*log10(abs(som_psd)) - maxdB1);
title('Densidade Espectral de Potência - Mensagem de Voz - dB');
xlabel('frequência (Hz)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reamostrar

new_sampling_rate = fs * os / downsampling_rate;

sampling_fun = [ 1; zeros( downsampling_rate - 1, 1 ) ];
sampling_fun = repmat( sampling_fun, length(som)/length(sampling_fun), 1 );
som_ds = som .* sampling_fun;

figure(3)
plot( t1, som_os, t1, som_ds );

% estimar densidade espectral de potência
[ som_ds_psd, f] = pwelch( som_ds, [], [],[], fs*os);

figure(4);
plot( f, abs(som_ds_psd) );
titlestr = sprintf('Densidade Espectral de Potência - Mensagem de Voz, fa = %f Hz', new_sampling_rate);
title(titlestr);
xlabel('frequência (Hz)');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/q
% filtragem



% definir filtro passa baixa
ffil = 4000; % frequência de corte
lfil = 100; % tamanho da resposta impulsional
h= fir1 ( lfil, ffil / (fs*os) * 2 );

% figure( 5 )
% freqz( h );

% filtrar
s_fil = filter( h, 1, som_ds * downsampling_rate );

energy1 = sum( som_ds .^2 );
energy_fil = sum( s_fil .^2 );
norm = sqrt( energy1 / energy_fil );

wavwrite(s_fil * norm , fs * os, bits, 'reamostrado.wav')


% estimar densidade espectral de potência
[ som1_f_psd ] = pwelch( s_fil, [], [], [], fs*os);

figure(5);
plot(f, abs(som_psd),f, abs(som1_f_psd));
title('Densidade Espectral de Potência - Mensagem de Voz');
xlabel('frequência (Hz)');
legend('original','filtrado');


