% demo_ampModulate - script to demonstrate ampModulate function
%
% Usage:
%   demo_ampModulate
%
% See also:
%   ampModulate

%% Demo 1: amplitude modulate a few cycles of a sinusoid
% choose frequencies that provide a nice plot
fm = 1;
fc = 8;
fs = 200;
BW = 1;

dur = 3/fm; % plot 3 cycles

% make a sinusoidal message signal
ttm = 0:1/fs:dur;
ss = 1*cos(2*pi*fm*ttm);

% amplitude modulation
[sp,ce] = ampModulate(ss, fc, BW, 0.5, fs);

% plot
figure(1)
ttp = (0:length(sp)-1)/fs; % sp is not the same length as ss due to filtering

plot(ttp,sp,'b', ttm, ss, 'r', ttp, ce, 'r--');
grid on
xlabel('Time (s)')
ylabel('Amplitude')
legend('Passband signal', 'Message Signal', 'Envelope', 'Location', 'SouthEast')

%% Demo 2: amplitude modulate a speech signal 
% load the speeach signal
[sss, fss] = wavread('kennedy_moon.wav');

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

% set bandwidth and carrier frequency
fcs = 8e3;
BWs = 2e3;

% amplitude modulation
[sps,ce] = ampModulate(sss, fcs, BWs, 0.5, fss);

% analyze the spectrum of the passband signal
h = spectrum.welch;
set(h, 'SegmentLength', 1024)
Hpsd = psd(h, sps, 'Fs', fss);         % Calculate the PSD 

figure(2)
plot(Hpsd)                           % Plot the PSD.

% annotate plot
hold on
yy = get(gca,'Ylim');
plot([1 1]*fcs/1e3, yy, 'b--', [1 1]*(fcs+BWs)/1e3, yy, 'r--', [1 1]*(fcs-BWs)/1e3, yy, 'r--')
hold off
legend('Signal', 'Carrier Freq', 'Carrier Freq \pm BW')
