% Neste script é mostrada a densidade espectral de potência de sinais de
% voz diferentes, e é mostrado o efeito de uma filtro passa-baixa nestes
% sinais.
%
% Universidade de Brasília
% André Noll Barreto, 11 set 2014


close all

%%%%%%%
% parâmetros
arquivoSom1 = 'female1.wav';
arquivoSom2 = 'male1.wav';
fmin = 32000; % frequência máxima observável na DFT
fc = 15000; %frequência da portadora 
ncorr = 256; % número de amostras da autocorrelação

%%%%%%%

[som1, fs, bits] = wavread(arquivoSom1);
[som2, fs, bits] = wavread(arquivoSom2);

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

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


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

t2 = (0:length(som2)-1) / fs;
t2_os = (0 : 1/os : (length(som2)-1)) / fs; % instantes de amostra
t2_s = 1 / fs / os;

som1_os = interp1(t1, som1, t1_os);  % som super-amostrado
som2_os = interp1(t2, som2, t2_os);  % som super-amostrado

% estimar densidade espectral de potência
[ som1_psd, f] = pwelch( som1_os, [], [],[], fs*os);
[ som2_psd ] = pwelch( som2_os, [], [], [], fs*os);

figure(1);
plot(f, abs(som1_psd),f, abs(som2_psd));
title('Densidade Espectral de Potência - Mensagem de Voz');
xlabel('frequência (Hz)');
legend('feminino','masculino');


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


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

figure(3);
plot(t1_os, som1_os);
title ('Sinal Mensagem');
xlabel('tempo (s)');


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

figure( 4 )
freqz( h );

% filtrar
s_fil = filter( h, 1, som1_os );

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

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


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

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


