% DSB_AM - demo script for conventional AM and envelope detection
 
%% AM parameters
fc = 5e3;    % carrier frequency
a_mod = 0.5; % modulation index

%% read in message signal, filter, and upsample it
[mm, fs] = audioread('kennedy_moon.wav');

% lowpass filter to limit to 2KHz
hh_mm = firpm(51, [0 1.5e3 2.0e3 fs/2]/(fs/2), [1, 1, 0, 0]);
mm = conv(mm,hh_mm);

% resample to higher sampling rate (4*fs)
mm = mm.';
mm = resample(mm, 4, 1);
fs = 4*fs;

% play sound and display spectrogram of message signal
spectrogram(mm,1024,0,1024,fs,'yaxis');
soundsc(mm,fs);

pause

%% DSB-AM
% carrier signal
tt = (0:length(mm)-1)/fs;
cc_t = cos(2*pi*fc*tt);

% DSB-AM signal with carrier amplitude set according to modulation index
Ac = abs(min(mm))/a_mod;
Ac = 0;
ss = (Ac + mm) .* cc_t;

% play sound and display spectrogram
spectrogram(ss,1024,0,1024,fs,'yaxis');
soundsc(ss,fs);

pause

plot(tt, ss, tt, (Ac + mm),'--',tt, -(Ac + mm),'--');
set(gca, 'XLim', 1 + [0, 30/fc]);
grid on
xlabel('Time (s)')
legend('Modulated Signal', 'Upper Envelope', 'Lower Envelope');

pause

%% Attempt Coherent Demodulation
% carrier at receiver with phase error
disp('Demod - no phase error')
theta = 0;                  % no phase error
% theta = pi/2;               % constant phase offset (worst case)
% theta = 2*pi*10*tt + pi/2;  % frequency offset of 10 Hz

cc = cos(2*pi*fc*tt + theta);

% multiply received signal with carrier
yy = ss.*cc;

disp('Before LPF')
spectrogram(yy,1024,0,1024,fs,'yaxis');
soundsc(yy,fs);

pause

% lowpass filter to remove double frequency components
hh_lpf = firpm(150, [0, 2e3, 3e3, fs/2]/(fs/2), [1, 1, 0, 0]);
mm_hat = conv(yy, hh_lpf);

disp('After LPF')
spectrogram(mm_hat,1024,0,1024,fs,'yaxis');
soundsc(mm_hat,fs);

pause

disp('Demod - worst case phase offset')
theta = pi/2;               % constant phase offset (worst case)
% theta = 2*pi*10*tt + pi/2;  % frequency offset of 10 Hz

cc = cos(2*pi*fc*tt + theta);

% multiply received signal with carrier
yy = ss.*cc;
mm_hat = conv(yy, hh_lpf);
spectrogram(mm_hat,1024,0,1024,fs,'yaxis');
soundsc(mm_hat,fs);

pause

disp('Demod - frequency offset')
% theta = pi/2;               % constant phase offset (worst case)
theta = 2*pi*0.5*tt + pi/2;  % frequency offset of 10 Hz

cc = cos(2*pi*fc*tt + theta);

% multiply received signal with carrier
yy = ss.*cc;
mm_hat = conv(yy, hh_lpf);
spectrogram(mm_hat,1024,0,1024,fs,'yaxis');
soundsc(mm_hat,fs);
