% FM.m
%
% Neste script são demonstrados os esquemas de modulação FM e PM
% É observado o espectro e o sinal no tempo para sinais modulante e
% modulado.
%
% A demonstração é feita com uma amostra de áudio gravada em arquivo tipo 
% wav sem compressão

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

%%%%%%%%%%%%%%
% parâmetros
arquivoSom = 'female2.wav';
fmin = 256000; % frequência máxima observável na DFT
fc = 15000; %frequência da portadora 
ncorr = 256; % número de amostras da autocorrelação
recuo = -3; % recuo de potência em dB
B = 4e3; % banda do sinal em Hz
filsize = 100;

SNR_dB = 0;
%%%%%%%%%%%%%%
SNR = 10^(SNR_dB/10);

[som, fs, bits] = wavread(arquivoSom);
pot_som = mean(abs(som).^ 2);
t = (0:length(som)-1).' / fs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% super amostrar sinal para aumentar faixa de frequência observável
os = ceil(2*fmin/fs); % fator de superamostragem

t_s = 1 / fs / os;
f_os = fs * os;

som_os = interp(som, os);  % som super-amostrado
t_os = (0:length(som_os)-1).' / f_os;


% estimar densidade espectral de potência
figure(1);
[PSD_som, f_som] = pwelch(som_os, [], [], [], f_os, 'twosided' );
PSD_som = 10*log10(PSD_som);
PSD_som = PSD_som - max(PSD_som);
plot(f_som / 1000, PSD_som);
xlabel('frequência(kHz)')
ylabel('dB/Hz')
title('Densidade Espectral de Potência - Mensagem de Voz');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% modulação AM DSB-SC
s_dsb =  som_os .* cos(2 * pi * fc * t_os);

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

[PSD_dsb, f_dsb] = pwelch(s_dsb, [], [], [], f_os, 'twosided' );
PSD_dsb = 10*log10(PSD_dsb);
PSD_dsb = PSD_dsb - max(PSD_dsb);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  não linearidade
clip_level = sqrt(mean(abs(s_dsb).^2) * 10^(recuo/10));
index = find(abs(s_dsb)>clip_level);
s_dsb_nl = s_dsb;
s_dsb_nl(index) = clip_level * sign(s_dsb(index));

% ruído
sigPower = mean( s_dsb .^2 );
noisePower = sigPower / SNR;
s_dsb_noise = s_dsb + randn( size( s_dsb ) ) * sqrt(noisePower);

figure(3);
plot(t_os, s_dsb, t_os, s_dsb_nl, t_os, s_dsb_noise);
legend('sinal AM-DSB original', 'sinal AM-DSB ceifado', 'sinal AM-DSB com ruído');

figure(4);
[PSD_dsb_nl, f_dsb_nl] = pwelch(s_dsb_nl, [], [], [], f_os, 'twosided' );
PSD_dsb_nl = 10*log10(PSD_dsb_nl);
PSD_dsb_nl = PSD_dsb_nl - max(PSD_dsb_nl);

[PSD_dsb_noise, f_dsb_noise] = pwelch(s_dsb_noise, [], [], [], f_os, 'twosided' );
PSD_dsb_noise = 10*log10(PSD_dsb_noise);
PSD_dsb_noise = PSD_dsb_noise - max(PSD_dsb_noise);


plot(f_dsb/1000, PSD_dsb, f_dsb_nl/1000, PSD_dsb_nl, f_dsb_noise/1000, PSD_dsb_noise)
title('Densidade Espectral de Potência');
xlabel('frequência (kHz)');
ylabel('dB');
legend('sinal AM-DSB original', 'sinal AM-DSB ceifado', 'sinal AM-DSB com ruído');




%%%%%%%%%%%%%%%%%%%%%%%%%%%
% demodulação
s_dsb_demod = 2* s_dsb .* cos(2 * pi * fc * t_os);
s_dsb_demod_nl = 2 * s_dsb_nl .* cos(2 * pi * fc * t_os);
s_dsb_demod_noise = 2 * s_dsb_noise .* cos(2 * pi * fc * t_os);


% definir filtro passa baixa
h = fir1 (filsize, B / f_os);

% filtrar
s_dsb_demod = filter(h,1,s_dsb_demod);
s_dsb_demod_nl = filter(h,1,s_dsb_demod_nl);
s_dsb_demod_noise = filter(h,1,s_dsb_demod_noise);

figure(5)
plot(t_os, som_os, t_os, s_dsb_demod, t_os, s_dsb_demod_nl, t_os, s_dsb_demod_noise  );
legend('sinal original', 'sinal demodulado', ...
       'sinal com distorção não-linear demodulado', ...
       'sinal com ruído');

aux = decimate(s_dsb_demod, os);
aux = aux * sqrt(pot_som / mean(abs(aux).^ 2));
wavwrite(aux, fs, bits, 'dsbsc_demod.wav')

aux = decimate(s_dsb_demod_nl, os);
aux = aux * sqrt(pot_som / mean(abs(aux).^2));
wavwrite(aux, fs, bits, 'dsbsc_demod_nl.wav')

aux = decimate(s_dsb_demod_noise, os);
aux = aux * sqrt(pot_som / mean(abs(aux).^2));
wavwrite(aux, fs, bits, 'dsbsc_demod_noise.wav')




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% modulação FM

a = cumsum(som_os) * t_s;
mp = max(abs(som_os));

kfmp = B/20;
kf = kfmp/ mp;


% banda estreita
s_fm_nb = cos(2*pi*fc*t_os) - kf * a.* sin(2*pi*fc*t_os);
[PSD_fm_nb f_fm_nb] = pwelch(s_dsb_nl, [], [], [], f_os, 'twosided' );
PSD_fm_nb = 10*log10(PSD_fm_nb);
PSD_fm_nb = PSD_fm_nb - max(PSD_fm_nb);

% banda larga
Bfm = 40e3;
kf = (Bfm / 2 - B) / mp;

s_fm_40 =  cos(2 * pi * fc * t_os + 2 * pi * kf * a);
s_fm40_noise = s_fm_40 + randn( size( s_fm_40 ) ) * sqrt(noisePower);


[PSD_fm_40 f_fm_40] = pwelch(s_fm_40, [], [], [], f_os, 'twosided' );
PSD_fm_40 = 10*log10(PSD_fm_40);
PSD_fm_40 = PSD_fm_40 - max(PSD_fm_40);

Bfm = 160e3;
kf = (Bfm / 2 - B) / mp;

s_fm_80 =  cos(2 * pi * fc * t_os + 2 * pi * kf * a);
s_fm80_noise = s_fm_40 + randn( size( s_fm_80 ) ) * sqrt(noisePower);

[PSD_fm_80 f_fm_80] = pwelch(s_fm_80, [], [], [], f_os, 'twosided' );
PSD_fm_80 = 10*log10(PSD_fm_80);
PSD_fm_80 = PSD_fm_80 - max(PSD_fm_80);

figure(6);
plot(t_os, s_fm_nb, t_os, s_fm_40);
xlabel('t');
legend( 'FM-NB', 'FM, kf = 5')

figure(7);
% plot( f_dsb/1000, PSD_dsb, f_fm_nb/1000, PSD_fm_nb, f_fm_40/1000, PSD_fm_40, ...
%       f_fm_80/1000, PSD_fm_80 );
  plot( f_dsb/1000, PSD_dsb,  f_fm_40/1000, PSD_fm_40, ...
      f_fm_80/1000, PSD_fm_80 );

title('Densidade Espectral de Potência - Sinal Modulado FM');
xlabel('frequência (kHz)');
ylabel('dB');
%legend('DSB-SC', 'FM banda estreita', 'FM, B_FM = 40kHz', 'FM, B_FM = 80kHz'); 
bfm_string =  sprintf( 'FM, B_FM = %dkHz',  Bfm/1000 );
legend('DSB-SC',  'FM, B_FM = 40kHz', bfm_string); 


%%%%%%%%%%%%%%%%%
% demodulação FM
% discriminador de frequência

s_fm_demod = diff(s_fm_80) / t_s;
s_fm_demod = [0;s_fm_demod];

% detector de envelope
s_fm_demod(s_fm_demod < 0) = 0; % diodo
s_fm_demod = filter(h,1,s_fm_demod); % filtro RC
s_fm_demod = s_fm_demod - mean(s_fm_demod);% bloqueador DC
s_fm_demod = s_fm_demod / max(abs(s_fm_demod)) / 2;

aux = decimate(s_fm_demod, os);
aux = aux * sqrt(pot_som / mean(abs(aux).^ 2));
wavwrite(aux, fs, bits, 'fm_demod.wav')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% discriminador de frequência 40 kHz + ruído
s_fm_demod = diff(s_fm_40) / t_s;
s_fm_demod = [0;s_fm_demod];

% detector de envelope
s_fm_demod(s_fm_demod < 0) = 0; % diodo
s_fm_demod = filter(h,1,s_fm_demod); % filtro RC
s_fm_demod = s_fm_demod - mean(s_fm_demod);% bloqueador DC
s_fm_demod = s_fm_demod / max(abs(s_fm_demod)) / 2;

aux = decimate(s_fm_demod, os);
aux = aux * sqrt(pot_som / mean(abs(aux).^ 2));
wavwrite(aux, fs, bits, 'fm_demod_40_noise.wav')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% discriminador de frequência 160 kHz + ruído
s_fm_demod = diff(s_fm_80) / t_s;
s_fm_demod = [0;s_fm_demod];

% detector de envelope
s_fm_demod(s_fm_demod < 0) = 0; % diodo
s_fm_demod = filter(h,1,s_fm_demod); % filtro RC
s_fm_demod = s_fm_demod - mean(s_fm_demod);% bloqueador DC
s_fm_demod = s_fm_demod / max(abs(s_fm_demod)) / 2;

aux = decimate(s_fm_demod, os);
aux = aux * sqrt(pot_som / mean(abs(aux).^ 2));
wavwrite(aux, fs, bits, 'fm_demod_160_noise.wav')


%%%%%%%%%%%%%%%%%%%%
% não linearidade
clip_level = sqrt(mean(abs(s_fm_80).^2) * 10^(recuo/10));
index = find(abs(s_fm_80)>clip_level);
s_fm_nl = s_fm_80;
s_fm_nl(index) = clip_level * sign(s_fm_80(index));

figure(8);
plot(t_os, s_fm_80, t_os, s_fm_nl);
legend('sinal FM original', 'sinal FM ceifado');

%%%%%%%%%%%%%%%%%
% demodulação FM
% discriminador de frequência

s_fm_nl_demod = diff(s_fm_nl) / t_s;
s_fm_nl_demod = [0;s_fm_nl_demod];

% detector de envelope
s_fm_nl_demod(s_fm_nl_demod < 0) = 0; % diodo
s_fm_nl_demod = filter(h,1,s_fm_nl_demod); % filtro RC
s_fm_nl_demod = s_fm_nl_demod - mean(s_fm_nl_demod);% bloqueador DC
s_fm_nl_demod = s_fm_nl_demod / max(abs(s_fm_nl_demod)) / 2;

aux = decimate(s_fm_nl_demod, os);
aux = aux * sqrt(pot_som / mean(abs(aux).^ 2));
wavwrite(aux, fs, bits, 'fm_nl_demod.wav')


% demodulação FM com ruído


