% demo_pulseShape - demonstrate the impact of pulse shape in time and
% frequency domain

T = 1;
fs = 16;
fsT = fs*T;

sigma=0.35;  % controls width of Gaussian pulse
aa = 0.5;    % excess bandwidth

%% menu
fprintf('\n\n Select one of the following pulses:\n');
fprintf('\t (1) Rectangular\n');
fprintf('\t (2) Half-sine\n');
fprintf('\t (3) Squared Half-sine\n');
fprintf('\t (4) Gaussian (sigma = %4.2f)\n', sigma);
fprintf('\t (5) Sinc\n');
fprintf('\t (6) Trapezoidal (excess BW: %d percent)\n', aa*100);
fprintf('\t (7) Raised Cosine (excess BW: %d percent)\n\n', aa*100);

choice = input('Choose 1-7: ');

%% pulse shape
switch choice
    % full response pulses
    case 1
        pulse = ones(1,fsT);             % rectangular
    case 2
        pulse = sin(pi*(0:fsT-1)/fsT);   % half sine
    case 3
        pulse = sin(pi*(0:fsT-1)/fsT).^2;% squared half sine

        % partial response pulses
    case 4
        ttp = -10:1/fsT:10-1/fsT;
        pulse = exp(-0.5/sigma^2*ttp.^2);        % Gaussian pulse - not Nyquist

        % Nyquist pulses
    case 5
        pulse = sinc(ttp);   % sinc pulse
    case 6
        pulse = sinc(ttp).*sinc(aa*ttp); % trapezoidal
    case 7
        pulse = rcos(ttp,aa);
end

figure(3)
plot((0:length(pulse)-1)/fsT, pulse)
grid on
xlabel('Time symbol periods')
ylabel('Amplitude');

%% make a short sequence of BPSK symbols
b = [1, -1, 1, 1, 1, 1, 1, -1, -1, -1];

ss = linearBBMod( b, pulse, fsT );
tt = (0:length(ss)-1)/fsT;

figure(1)
sampTimes = length(pulse)/2+1 + (0:fsT:fsT*length(b)-1);
plot(tt, ss, 'b', ...
    tt(sampTimes), ss(sampTimes), 'ro',...
    tt(sampTimes), b, 'bx',...
    tt(sampTimes+1), ss(sampTimes+1), 'gd')
grid on
xlabel('Time (symbol periods)')
ylabel('Amplitude') 

legend('TX Signal', 'Samples at mT', 'Symbols at mT', 'Samples at mT+\delta')


%% part 2: estimate the PSD using a long symbol sequence
Nsymbols = 8196;
win = hamming(1024).';

pulse = pulse/norm(pulse);

% random symbols
symbols = 2*(rand(1,Nsymbols)>0.5) - 1;  % BPSK

% digital modulation
sig = linearBBMod(symbols, pulse, fsT);

[Sx, ff] = estimatePSD(sig, win, fs);

figure(2)
plot(ff, 10*log10(Sx));
grid on
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

%% Quantitative measurements
% ISI measures
sampTimes = length(pulse)/2+1 + (0:fsT:fsT*length(symbols)-1);
ISI = max(abs(sig(sampTimes) - symbols));
ISI_delta = max(abs(sig(sampTimes+1) - symbols));

% 99% power containment
NFFT = 2048;
gamma = 0.99;
PP = abs(fftshift(fft(pulse, NFFT)));
ff = (-NFFT/2:NFFT/2-1)*fsT/NFFT;

SSr = cumsum(PP.^2)/sum(PP.^2);
indB = find(SSr >= (1+gamma)/2, 1);
Bgamma = 2*ff(indB);

fprintf('Maximum ISI: %6.4f\n', ISI);
fprintf('Maximum ISI with timing error: %6.4f\n', ISI_delta);
fprintf('%d-percent containment BW: %6.4f\n\n', 100*gamma, Bgamma);
