% Signal space and constellation
% Ex 5
clc; clear all; 
% close all

% Input papameters
% For passband implemetation sampling frequency is a necessary parameter
fs = 1000; % sampling frequency [Hz]
Tsamp = 1/fs; % sampling time [s]
rb = 100; % bit rate [bit/sec]
N = 100; % number of bits to transmit
pulse_flag = 2; % 0 = rect, 1 = rc, 2 = rrc

% Constellation or bit to symbol mapping
constellation = [(1 + 1i) (1 - 1i) (-1 -1i) (-1 + 1i)]/sqrt(2); % Constellation 1 - QPSK/4-QAM

%s = exp(1i*((0:3)*pi/2 + pi/4)); % Constellation 1 - same constellation generated as PSK
scatterplot(constellation); grid on;                            % Constellation visualization

M = length(constellation);                                      % Number of symbols in the constellation
bits_per_symb = log2(M);                                        % Number of bits per symbol
fd = rb/bits_per_symb;                                          % Symbol rate [symb/s]
fsfd = fs/fd;                                       % Number of samples per symbol (choose fs such that fsfd is an integer for simplicity) [samples/symb]
a = randsrc(1,N,[0 1]);                             % Information bits
m = buffer(a, bits_per_symb)';                      % Group bits into messages consisting of bits
sym_idx = bi2de(m, 'left-msb')'+1;                   % Bits to symbol index
x = constellation(sym_idx);                           % Look up symbols using the indices

x_upsample = upsample(x, fsfd);                     % Space the symbols fsfd apart, to enable pulse shaping using conv.

% Now we have a signal ready to be shaped.

switch pulse_flag
    case 0 %Rectangular
        pulse = ones(1, fsfd);                   % Generate a rectangular pulse
        pulse_train = conv(pulse,x_upsample);    % Each symbol replaced by the pulse shape and added
        
    case 1 %RC
        span = 6;
        [pulse, t] = rcpuls(0.35,1/fd,fs,span);
        pulse_train = conv(pulse,x_upsample);
        
    case 2 %RRC
        span = 6;
        [pulse, t] = rtrcpuls(0.4,1/fd,fs,span);
        pulse_train = conv(pulse,x_upsample);
        
    otherwise
        disp('Unknown pulse');
        return;
end

figure; subplot(2,1,1); plot( real(pulse_train), 'b');
title('real')
subplot(2,1,2); plot(Tsamp*(0:(length(pulse_train)-1)), imag(pulse_train), 'b');
title('imag')

% END OF EXERCISE 5
%% EX. 6

% Add AWGN noise to the signal
SNRdB = 20; %decide noise level
y = awgn(pulse_train, SNRdB, 'measured'); % add noise

% Plot the received signal 'y'
figure; subplot(2,1,1); plot(Tsamp*(0:(length(y)-1)), real(y), 'b');
title('real')
subplot(2,1,2); plot(Tsamp*(0:(length(y)-1)), imag(y), 'b');
title('imag')


MF = fliplr(conj(pulse));        %create matched filter impulse response
MF_output = filter(MF,1,y);      % run received signal through matched filter
figure; plot(real(MF_output))
MF_output = MF_output(length(MF):end); %remove transient
MF_output_conv = conv(pulse, y);  % Another approach to MF using conv, what's the difference?
figure; plot(real(MF_output))
MF_output_conv = MF_output_conv(length(MF):end-length(MF)+1);

rx_vec = MF_output(1:fsfd:end);  %get sample points (assuming perfect synchronisation).

%%%%%%%%% MY STUFF %%%%%%%%%%%%%
% %% Introduce phase offset in received baseband signal.
% MF_output_offset = circshift(MF_output, 5);
% rx_vec = MF_output_offset(1:fsfd:end);      % Sample at given phase offset.
% 
% sync_offset = 0; % Samples to offset by in order to have synchronization.
% %figure; plot(real(MF_output_conv))
% 
% scatterplot(rx_vec); %scatterplot of received symbols
% 
% figure
% hold off
% avg = zeros(fsfd, 1);
% for delta = 1:fsfd
%     shifted_MF_output = circshift(MF_output_offset, delta-1);
%     shifted_samples = shifted_MF_output(1:fsfd:end);
%     avg(delta) = sum(abs(shifted_samples))./length(shifted_samples);
%     
%     plot(shifted_samples, 'ko')
%     xlabel('In Phase'); ylabel('Quadrature');
%     title(['Sync offset: ',  num2str(delta-1)])
%     %pause(0.5)
% end
% 
% [max_avg, max_idx] = max(avg);
% sync_offset = max_idx-1;
% 
% figure
% plot(0:fsfd-1,avg)
% xlabel('sync offset [samples]'); ylabel('Average energy');
% 
% MF_output_corrected = circshift(MF_output_offset, sync_offset);
% scatterplot(MF_output_corrected(1:fsfd:end))
% title(['sync offset ', num2str(sync_offset)]) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure
subplot(2,1,1)
plot(real(y(span*fsfd:end-span*fsfd)));
subplot(2,1,2);
plot(real(MF_output));
hold on; stem((1:fsfd:length(MF_output)),real(MF_output(1:fsfd:end)));

% Minimum Eucledian distance detector
% Relate the detection to Detection region
metric = abs(repmat(rx_vec.',1,4) - repmat(constellation, length(rx_vec), 1)).^2; % compute the distance to each possible symbol
[tmp m_hat] = min(metric, [], 2); % find the closest for each received symbol
m_hat = m_hat';   % get the index of the symbol in the constellation

SER = symerr(sym_idx, m_hat) %count symbol errors
b_hat_buffer = de2bi(m_hat-1, 2, 'left-msb')'; %make symbols into bits
b_hat = b_hat_buffer(:)'; %write as a vector
BER = biterr(a, b_hat) %count of bit errors



%% EX. 7 
eyediagram(pulse_train(span*fsfd:end-span*fsfd), fsfd, 1/fd) % plot the eyediagram from the output of pulse shaping
eyediagram(MF_output, fsfd, 1/fd); % plot the eyediagram from the output of matched filter using MATLAB's function
eyediagram(pulse_train(span*fsfd+5:end-span*fsfd), fsfd, 1/fd) % sychronization error, i.e., wrong sampling time
eyediagram(MF_output.*exp(1j*0.4), fsfd, 1/fd); %phase error
eyediagram(abs(MF_output), fsfd, 1/fd);
eyediagram(abs(MF_output.*exp(1j*0.4)), fsfd, 1/fd);

% write your own function to plot eyediagram
MF_eye = reshape(MF_output, fsfd, N/bits_per_symb); %reshape in terms of samples
figure; plot(circshift(real(MF_eye), fsfd/2)); 
figure;plot(real(buffer(MF_output((span-0.5)*fsfd+1:end-(span+0.5)*fsfd), fsfd))) %plot eye diagram, shift fsfd/2 to find sampling instant
freqz(pulse_train)

% Plot frequency response of the pulse train, is it bandlimited and has a BW corresponding to the pulse chosen?
% Write your own function for plotting eye diagram
% What does a phase shift do to the eyediagram?
% Phase has to be determined in a receiver (part of synchronization)
% Eye diagram of abs(pulse_tr) in case of phase shift, is this useful?

