% QAM.m
%
% Neste script é demonstrado o esquema de modulação QAM, onde modulamos 
% dois sinais diferentes em uma mesma banda passante.
% São observados os espectros dos sinais modulante e modulado.
% Também é demodulado o sinal, e mostrado que para sistemas DSB-SC desvios de
% fase podem causar interfeência entre os canais I e Q.
%
% A demonstração é feita com amostras 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:
%   "testqam1_demod.wav"
%   "testqam2_demod.wav" - sinais demodulados perfeitamente
%   "testqam1_demod_err.wav"
%   "testqam2_demod_err.wav" - sinais demodulados com erro de fase

% Universidade de Brasília
% André Noll Barreto, 16 jun 2010

close all % fechar janelas abertas

%%%%%%%%%%%%%%%%
% parâmetros
arquivoSom = 'female2.wav';
arquivoSom1 = 'male1.wav';

fmin = 32000; % frequência máxima observável na DFT
fc = 15000; %frequência da portadora 
phi_err = pi / 8; % erro de fase na demodulação DSB-SC
ncorr = 1024; % número de amostras da FFT
freq_err = 0.01; % desvio de frequência em relação a fc
%%%%%%%%%%%


[som, fs, bits] = wavread(arquivoSom);
[som1, fs1, bits1] = wavread(arquivoSom1);

if (fs1 ~= fs | bits1 ~= bits)
  errormsg('amostras de áudio devem ter mesmos parâmetros');
end

minsize = min(length(som), length(som1));
som = som(1:minsize);
som1 = som1(1:minsize);

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

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

f0 = 1 / (ncorr * 2 * t_s);
f = (0:(ncorr*2-1)) * f0;

som_os = interp1(t, som, t_os);
som_os1 = interp1(t, som1, t_os);

[ som_psd, f ] = pwelch(som_os, ncorr, [],[], fs*os);
[ som_psd1 ] = pwelch(som_os1, ncorr, [],[], fs*os);

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

clear som som1 som_psd som_psd1;

s_qam =  som_os .* cos(2 * pi * fc * t_os) + som_os1 .* sin(2 * pi * fc * t_os);
qam_psd = pwelch (s_qam, ncorr, [],[], fs*os);

figure(2);
plot(f, qam_psd);
title('Densidade Espectral de Potência - Mensagem Modulada');
xlabel('frequência (Hz)');

clear qam_psd som_os som_os1;

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

% filtrar
som_demod = 2 * filter( h, 1, s_qam .* cos(2 * pi * fc * t_os) );

wavwrite(som_demod, fs * os, bits, 'testqam1_demod.wav');

som1_demod = 2 * filter( h, 1, s_qam .* sin(2 * pi * fc * t_os) );

wavwrite(som1_demod, fs * os, bits, 'testqam2_demod.wav');


% demodulação com erro de fase
theta = phi_err;

som_demod = 2 * filter( h, 1, s_qam .* cos(2 * pi * fc * t_os + theta) );
wavwrite(som_demod, fs * os, bits, 'testqam1_demod_err_theta.wav');

som1_demod = 2 * filter( h, 1, s_qam .* sin(2 * pi * fc * t_os + theta) );
wavwrite(som1_demod, fs * os, bits, 'testqam2_demod_err_theta.wav');

% demodulação com erro de frequência
som_demod = 2 * filter( h, 1, s_qam .* cos(2 * pi * fc*(1+freq_err) * t_os + theta) );
wavwrite(som_demod, fs * os, bits, 'testqam1_demod_err_freq.wav');

som1_demod = 2 * filter( h, 1, s_qam .* sin(2 * pi * fc* (1+freq_err) * t_os + theta) );
wavwrite(som1_demod, fs * os, bits, 'testqam2_demod_err_freq.wav');


