%
% Neste script são demonstrados os esquemas de modulação AM (DSB+C) e DSB-SC.
% É observado o espectro e o sinal no tempo para sinais modulante e modulado.
% Também é demodulado o sinal, e mostrado que para sistemas DSB-SC desvios de
% fase e frequência podem atenuar e /ou distorcer o sinal.
%
% A demonstração é feita com uma amostra de áudio gravada em arquivo tipo wav
% obs: arquivos wav com compressão não podem ser abertos pelo SciLab
%
% são salvos os arquivos de áudio:
%   "dsbsc_demod.wav" - sinal DSB-SC demodulado perfeitamente
%   "dsbsc_demod_err_phi.wav" - sinal DSB-SC demodulado com erro de fase
%   "dsbsc_demod_err_freq.wav" - sinal DSB-SC demodulado com erro de fase
%   "limitado.wav" - sinal de áudio limitado em 50% da amplitude máxima

% Universidade de Brasília
% André Noll Barreto, 07 mai 2015 - convertido para Matlab

close all  % fechar janelas abertas

%%
% parâmetros
arquivoSom = 'female2.wav';
fmin = 32000; % frequência máxima observável na DFT
fc = 15000; % frequência da portadora 
ncorr = 1024; % número de amostras da FFT
phi_err = pi / 4; % erro de fase na demodulação DSB-SC
freq_err = 0.01; % erro de frequÃ¢ncia (relativo Ã  portadora)


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

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

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

t = (0:length(som)-1) / fs;
t_os = (0 : 1/os : (length(som)-1)) / fs; % instantes de amostra
t_s = 1 / fs / os;

som_os = interp1(t, som, t_os);  % som super-amostrado

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


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


subplot(1,2,2)
plot(f, 10*log10(abs(som_psd)));
title('Densidade Espectral de Potência (dB) - Mensagem de Voz');
xlabel('frequência (Hz)');
v = axis;
v(2) = 20000;
axis(v);


%%
% modulação AM DSB-SC

s_dsb_sc =  som_os .* cos(2 * pi * fc * t_os);

figure(2);
plot(t_os, som_os, t_os, s_dsb_sc);
title ('Sinal Mensagem e Sinal Modulado DSB-SC');
xlabel('tempo (s)');
legend('sinal mensagem', 'sinal modulado');

figure(3);
mod_psd = pwelch( s_dsb_sc, ncorr, [],[], fs*os);
plot(f, 10*log10(abs(mod_psd)));
title('Densidade Espectral de Potência - Sinal Modulado DSB-SC');
xlabel('frequência (Hz)');

% %
% demodulação perfeita

% multiplicar por cos( 2 pi fc t)
s_dsb_sc_demod = s_dsb_sc .* cos(2 * pi * fc * t_os);

figure(4);
demod1_psd = pwelch( s_dsb_sc_demod, ncorr, [],[], fs*os); 
plot(f, abs(demod1_psd));
title('Densidade Espectral de Potência - Sinal Demodulado DSB-SC (pré filtragem)');
xlabel('frequência (Hz)');

% definir filtro passa baixa
ffil = fs / 2;
lfil = 100;
h= fir1 ( lfil, ffil / (fs*os) * 2 );
 
% filtrar
s_dsb_sc_demod = filter(h, 1, s_dsb_sc_demod );

figure(5);
demod2_psd = pwelch( s_dsb_sc_demod, ncorr, [],[], fs*os);
plot(f, abs(demod2_psd));
title('Densidade Espectral de Potência - Sinal Demodulado DSB-SC (pós filtragem)');
xlabel('frequência (Hz)');

wavwrite(s_dsb_sc_demod, fs * os, bits, 'dsbsc_demod.wav')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% demodulaÃ§Ã£o com erro de fase

s_dsb_sc_demod = s_dsb_sc .* cos(2 * pi * fc * t_os + phi_err);
s_dsb_sc_demod = filter( h, 1, s_dsb_sc_demod );

figure(6);
plot(t_os, som_os, t_os, s_dsb_sc_demod );
title ('Sinal Mensagem Original e Sinal DeModulado DSB-SC com Erro de Fase');
xlabel('tempo (s)');
legend('sinal mensagem', 'sinal demodulado com erro de fase');

wavwrite(s_dsb_sc_demod, fs * os, bits, 'dsbsc_demod_err_phi.wav')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% demodulação com erro de frequência
 

 
s_dsb_sc_demod = s_dsb_sc .* cos(2 * pi * fc * (1+freq_err) * t_os);
 
% filtro passa baixa
s_dsb_sc_demod = filter(h, 1, s_dsb_sc_demod );
figure(7);
plot(t_os, som_os, t_os, s_dsb_sc_demod );
title ('Sinal Mensagem Original e Sinal DeModulado DSB-SC com Erro de Frequência');
xlabel('tempo (s)');
legend('sinal mensagem', 'sinal demodulado com erro de frequência');
 
wavwrite(s_dsb_sc_demod, fs * os, bits, 'dsbsc_demod_err_freq.wav');
 
figure(8)
demod_errf_psd = pwelch( s_dsb_sc_demod, ncorr, [],[], fs*os);
plot(f, abs(som_psd), f, abs(demod_errf_psd));
title('Densidade Espectral de Potência - Sinal Demodulado DSB-SC (pós filtragem)');
xlabel('frequência (Hz)');

%%
% modulação AM (DSB + C)

mmax = max(abs(som_os));
s_am =  2*(mmax + som_os) .* cos(2 * pi * fc * t_os);

eficiencia = mean(som_os.^2) / mean(s_am.^2);

titlestr = ['Sinal Modulado AM (eficiência = ' num2str(eficiencia) ')'];
figure(9);
plot(t_os, som_os, t_os, s_am);
title (titlestr);
xlabel('tempo (s)');
legend('sinal mensagem', 'sinal modulado');

figure(10);
mod_psd = pwelch( s_am, ncorr, [],[], fs*os);
plot(f, 10*log10(abs(mod_psd)));
title('Densidade Espectral de Potência - Sinal Modulado AM');
xlabel('frequência (Hz)');

% modulação com limitador de amplitude
lim = 0.1;
indlim = find(abs(som_os) > lim*mmax);
som_os_lim = som_os;
som_os_lim(indlim) = lim * mmax .* sign(som_os(indlim));
wavwrite(som_os_lim, fs * os, bits, 'limitado.wav');

mmax = max(abs(som_os_lim));
s_am =  2*(mmax + som_os_lim) .* cos(2 * pi * fc * t_os);

eficiencia = mean(som_os_lim.^2) / mean(s_am.^2);

titlestr = ['Sinal Modulado AM (eficiência = ' num2str(eficiencia) ') - com limitador'];
figure(11);
plot(t_os, som_os_lim, t_os, s_am);
title (titlestr);
xlabel('tempo (s)');
legend('sinal mensagem', 'sinal modulado');

figure(12);
mod_psd =  pwelch( s_am, ncorr, [],[], fs*os);
plot(f, 10*log10(abs(mod_psd)));
title('Densidade Espectral de Potência - Sinal Modulado AM - com limitador');
xlabel('frequência (Hz)');
