% script for HW problem 4.17

%% parameters
fsT = 32;
T = 1;
fs = fsT/T;

% half since pulse
hs_pulse = sqrt(2)*sin(pi*(0:fsT-1)/fsT);

%% part a: plot a short OQPSK signal
N = 10;

% create BPSK symbols for I and Q channels
bI = 2*(rand(1,N) > 0.5) - 1;
bQ = 2*(rand(1,N) > 0.5) - 1;

% linear baseband modulation
sI = linearBBMod(bI, hs_pulse, fsT);
sQ = linearBBMod(bQ, hs_pulse, fsT);

% turn into a complex signal; shift the Q signal by T/2
ss = complex(sI, [zeros(1,fsT/2), sQ(1:end-fsT/2)]);

% plot real and imaginary parts
figure(1)
subplot(2,1,1)
plot((0:N*fsT-1)/fsT, real(ss));
grid on
ylabel('In-phase signal');
title('I and Q signals are offset by T/2')

subplot(2,1,2)
plot((0:N*fsT-1)/fsT, imag(ss));
grid on
ylabel('Quadrature signal');
xlabel('Time (t/T)');

%% part b: plot magnitude and phase of complex baseband signal
figure(2)
subplot(2,1,1)
plot((0:N*fsT-1)/fsT, abs(ss));
grid on
ylabel('Magnitude');
title('Signal has constant Magnitude and linear phase')

subplot(2,1,2)
plot((0:N*fsT-1)/fsT, angle(ss)*180/pi);
grid on
ylabel('Phase (degrees)');
xlabel('Time (t/T)');

%% part c: plot the PSD of MSK signal
% PSD equals 2/T*|P(f)|^2 with |P(f)|^2 given by (4.13) in text
df = 0.001;
ff = -10:df:10;
SS = 2/T*8/pi^2*cos(pi*ff).^2./(1-4*ff.^2).^2;

% fix "gaps"
nanInd = find(abs(abs(ff)-0.5) < sqrt(eps));
SS(nanInd) = 2/T*0.5;

figure(3)
plot(ff, 10*log10(SS));
grid on
xlabel('Normalized Frequency (fT)');
ylabel('Power Spectral density (dB/Hz)');
axis([ff(1), ff(end), -80, 10]);

%% part e: compare MSK and OQPSK with rectangular pulses
SSr = 2/T*sinc(ff).^2;

plot(ff, 10*log10(SS/(df*sum(SS))), ff, 10*log10(SSr/(df*sum(SSr))));
grid on
xlabel('Normalized Frequency (fT)');
ylabel('Power Spectral density (dB/Hz)');
axis([ff(1), ff(end), -80, 10]);
legend('MSK', 'OQPSK - rect.')

%% part d): bandwidths
SSc = cumsum(SS)/sum(SS);
indBW = find(SSc >= (1+0.99)/2,1);
BW99 = 2*ff(indBW);

SSrc = cumsum(SSr)/sum(SSr);
indBW = find(SSrc >= (1+0.99)/2,1);
BW99r = 2*ff(indBW); % this value is not accurate unless max(ff) is VERY large

fprintf('MSK - 99%% containment BW:           %5.2f * 1/T\n', BW99);
fprintf('Minimum Nyquist BW:                 %5.2f * 1/T\n', 1);
fprintf('OQPSK (rect.) - 99%% containment BW: %5.2f * 1/T\n', 2*10.2);
