% Example of digital modulation and demodulation with BPSK.
%
% Rich Kozick, ELEC 470, Spring 1998

dt = 1e-4;           % Sample spacing (seconds)
Fs = 1/dt;           % Sampling rate
T = 1;               % Duration of signal to generate
t = (dt:dt:T)';      % Sample times
Nt = length(t);      % Number of sample points

Tb = 0.01;           % Bit period
Nb = T/Tb;           % Number of bits

SNR = 10;           % Signal to noise ratio (not dB)

% Generate bits

b = rand(Nb,1);
z0 = find(b < 0.5);
z1 = find(b >= 0.5);
b(z0) = -1*ones(size(z0));
b(z1) = +1*ones(size(z1));
m = kron(b,ones(Nt/Nb,1));

% Define carrier frequency and amplitude for modulation

fc = 1000;                  % Carrier frequency, in hertz
Ac = 1;
c = Ac * cos(2*pi*fc*t);    % Carrier wave, unmodulated
q = Ac * sin(2*pi*fc*t);

% Modulate the message, with added noise

An = Ac/sqrt(2*SNR);
nI = An*randn(Nt, 1);
nQ = An*randn(Nt, 1);
s = (m+nI) .* c - (nQ).*q;

figure(1)
plot(t,s);
title('MODULATED SIGNAL s(t)')
xlabel('Time')
ylabel('Amplitude')
axis([0 T -1.5*Ac 1.5*Ac]);

% Coherent demodulator

phi = 0;     % Phase of local oscillator
cI = 2*cos(2*pi*fc*t + phi);
cQ = 2*sin(2*pi*fc*t + phi);

vI = s .* cI;
vQ = s .* cQ;

% 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
[bb,aa] = butter(order, Fdig);

mI = filter(bb, aa, vI);   % Apply the filter
mI = mI - mean(mI);
mQ = -filter(bb, aa, vQ);   % Apply the filter
mQ = mQ - mean(mQ);

figure(2)
subplot(211)
plot(t,mI)
title('COHERENT DETECTOR OUTPUT: In-phase')
xlabel('Time')
subplot(212)
plot(t,mQ)
title('COHERENT DETECTOR OUTPUT: Quadrature')
xlabel('Time')

% Sample the demodulated signal to make decisions

NTb = Tb/dt;        % Number of samples per bit period
xIT = mI(NTb-1:NTb:Nt);
xQT = mQ(NTb-1:NTb:Nt);

% Create scatter plot of received constellation in the I-Q plane

figure(3)
plot(xIT, xQT, 'o')
axis([-1.5*Ac 1.5*Ac -1.5*Ac 1.5*Ac])
title('RECEIVED CONSTELLATION, *NO* MATCHED FILTER')
xlabel('In-Phase')
ylabel('Quadrature')

for k=1:Nb
  xITmf(k) = sum(mI( (k-1)*NTb+1:k*NTb ))/NTb;
  xQTmf(k) = sum(mQ( (k-1)*NTb+1:k*NTb ))/NTb;
end

figure(4)
plot(xITmf, xQTmf, 'o')
axis([-1.5*Ac 1.5*Ac -1.5*Ac 1.5*Ac])
title('RECEIVED CONSTELLATION WITH MATCHED FILTER')
xlabel('In-Phase')
ylabel('Quadrature')
