% Examples of modulation and demodulation of signals in MATLAB.
% Modulation types include AM, DSB-SC, and FM.
% The extension to SSB and VSB is simple.
% Use the imzoom command to zoom-in on interesting parts of the plots.
%
% Rich Kozick, ELEC 470, Spring 1998

dt = 1e-4;           % Sample spacing (seconds)
Fs = 1/dt;           % Sampling rate
t1 = (dt:dt:1)';    
t2 = (1+dt:dt:2)';
t = [t1; t2];        % Sample times: total of 2 seconds of data
Nt = length(t);      % Number of sample points
Nf = 2^18;           % Number of points in FFT, for faster computation
f = (0:Nf-1)'/Nf*Fs;   % Frequencies in FFT

% Construct message signal

rpp = (1:250)';      % Pieces of a triangle wave
rpn = (249:-1:0)';
rnn = -rpp;
rnp = -rpn;
tri = [rpp; rpn; rnn; rnp]/10;
tri10 = tri;
for k=2:10           % 10 cycles of triangle wave
  tri10 = [tri10; tri];
end

m = [tri10; 100*sinc(100*(t2-1.5))];    % Message
m = m - mean(m);                        % Remove DC
M = fft(m, Nf);
clear t1 t2 tri tri10 rnn rnp rpn rpp

figure(1)
subplot(211)
plot(t,m)
title('Message signal m(t)')
xlabel('Time (sec)')
subplot(212)
plot(f, abs(M))     
ax = axis;
axis([0 Fs/2 ax(3) ax(4)]);     % Plot up to one-half the sampling rate
title('Message spectrum M(f)')
xlabel('Frequency (Hz)')

clear M

% Define carrier frequency and amplitude for modulation

fc = 1000;
Ac = 1;
c = Ac * cos(2*pi*fc*t);    % Carrier wave, unmodulated

% AM modulation

ka = 0.009;

sAM = (1 + ka*m) .* c;
SAM = fft(sAM, Nf);

figure(2)
subplot(221)
plot(t,sAM)
title('AM modulated signal')
xlabel('Time (sec)')
subplot(223)
plot(f, abs(SAM))
ax = axis;
axis([0 Fs/2 ax(3) ax(4)]);
title('Spectrum of AM signal')
xlabel('Frequency (Hz)')

% Envelope detector for AM

sAMplus = hilbert(sAM);                  % Pre-envelope = s(t) + j * \hat{s}(t)
sAMtilde = sAMplus .* exp(-j*2*pi*fc*t); % Complex envelope
mAMdemod = abs(sAMtilde);                % Envelope detector output
mAMdemod = mAMdemod - mean(mAMdemod);    % Remove DC from output
MAMdemod = fft(mAMdemod, Nf);

subplot(222)
plot(t,mAMdemod)
title('Envelope detector output')
xlabel('Time (sec)')
subplot(224)
plot(f, abs(MAMdemod))
ax = axis;
axis([0 Fs/2 ax(3) ax(4)]);
title('Spectrum of envelope detector output')
xlabel('Frequency (Hz)')

clear sAM; clear SAM; clear SAMplus; clear SAMtilde;
clear mAMdemod; clear MAMdemod;

fprintf('QUESTION:  What if you reduce ka in the AM signal?\n\n');

% DSB-SC modulation

sDSB = m .* c;
SDSB = fft(sDSB, Nf);

figure(3)
subplot(221)
plot(t,sDSB)
title('DSB-SC modulated signal')
xlabel('Time (sec)')
subplot(223)
plot(f, abs(SDSB))
ax = axis;
axis([0 Fs/2 ax(3) ax(4)]);
title('Spectrum of DSB-SC signal')
xlabel('Frequency (Hz)')

% Coherent demodulator

phi = pi/6;     % Phase of local oscillator
local = cos(2*pi*fc*t + phi);

vDSB = sDSB .* local;

% Design a low-pass filter to remove components near 2*fc Hz.
% Use a Butterworth digital filter:
%    fcut = -3 dB cutoff frequency in Hz
%    order = filter order, which can be quite high in a digital filter.
%    Higher order means better filtering, but more computation.
% You can view the frequency response with:
%    freqz(b,a)
% In the plot from freqz, the frequency 1 corresponds to Fs/2 Hz = 5000 Hz.

Fcut = 1000;
order = 5;

Fdig = Fcut / (Fs/2);   % "Digital frequency" normalized by sampling rate
[b,a] = butter(order, Fdig);
% freqz(b,a)
% title('Filter for DSB-SC coherent detector')

% Note that you can design a high-pass filter similarly:
%   [b,a] = butter(order, Fdig, 'high');
% Then use the filter command to apply the filter to a signal.

mDSBdemod = filter(b, a, vDSB);   % Apply the filter to the vDSB signal
mDSBdemod = mDSBdemod - mean(mDSBdemod);
MDSBdemod = fft(mDSBdemod, Nf);

subplot(222)
plot(t,mDSBdemod)
title('Coherent detector output')
xlabel('Time (sec)')
subplot(224)
plot(f, abs(MDSBdemod))
ax = axis;
axis([0 Fs/2 ax(3) ax(4)]);
title('Spectrum of coherent detector output')
xlabel('Frequency (Hz)')

clear sDSB SDSB local vDSB mDSBdemod MDSBdemod

fprintf('QUESTION:  What happens when you change phi in the DSB-SC coherent detector?\n\n')
fprintf('QUESTION:  What is the best value for phi in this DSB-SC coherent detector?\n\n')
fprintf('QUESTION:  How would you demodulate SSB?\n\n')

% Frequency Modulation (FM)

kf = 2;      % Frequency sensitivity

int_m(1) = 0;    % Integrate the signal to produce FM signal
for k=1:Nt-1
  int_m(k+1) = int_m(k) + m(k)*dt;
end

sFM = Ac * cos(2*pi*fc*t + 2*pi*kf*int_m');
SFM = fft(sFM, Nf);

figure(4)
subplot(221)
plot(t,sFM)
title('FM modulated signal')
xlabel('Time (sec)')
subplot(223)
plot(f, abs(SFM))
ax = axis;
axis([0 Fs/2 ax(3) ax(4)]);
title('Spectrum of FM signal')
xlabel('Frequency (Hz)')

% FM demodulator:  Derivative followed by envelope detector

dsFM = sFM(1);   % Compute derivative
for k=2:Nt
  dsFM(k) = (sFM(k)-sFM(k-1))/dt;
end

% Envelope detector

dsFMplus = hilbert(dsFM);                      % Pre-envelope
dsFMtilde = dsFMplus .* exp(-j*2*pi*fc*(t'));  % Complex envelope
mFMdemod = abs(dsFMtilde);                     % Envelope detector output
mFMdemod = mFMdemod - mean(mFMdemod);          % Remove DC from output
MFMdemod = fft(mFMdemod, Nf);

subplot(222)
plot(t,mFMdemod)
title('FM detector output')
xlabel('Time (sec)')
subplot(224)
plot(f, abs(MFMdemod))
ax = axis;
axis([0 Fs/2 ax(3) ax(4)]);
title('Spectrum of FM detector output')
xlabel('Frequency (Hz)')

fprintf('QUESTION:  For FM, what is the frequency deviation?\n\n');
fprintf('QUESTION:  How does the FM signal bandwidth compare with Carsons rule?\n\n');
fprintf('QUESTION:  For FM, try changing kf, and see effects on the spectrum.\n\n');