%% Simulates GFSK communication and plots the performance in terms of BER
%
% For the CFO correction, the samples are processed "one by one" to match 
% the real behavior of the radio.
%

%% Transmission

%%% Transmission/modulation parameters

f_c = 2.45e09;          % Carrier frequency
r_symb = 1e06;          % Symbol rate
Ts = 1/r_symb;          % Symbol period
m = 0.5;                % Modulation index
BTb = 0.5;              % Bandwidth - Symbol perdiod product


N_symb = 1e06;          % Number of symbols used for the simulation
reso = 8;               % Resolution factor (number of samples by symbol period at TX)
dfppm = 30;				% CFO expressed in ppm wrt the carrier frequency


%%% For BER vs Eb/N0 plot and noise variance (N0) generation
Nsnr = 16;              % Number of different values for the plot
Eb_N0_dB = linspace(0,Nsnr-1,Nsnr);
Eb_N0 = 10.^(Eb_N0_dB/10);
Eb = 1;
N0 = Eb./Eb_N0;


%%% Generation of the preamble + noise before packet parameter
preamble = [1;-1;1;-1;1;-1;1;-1];
Lpreamble = length(preamble);

Lnoisestart = 92;		% Length of the noise portion preceding the packet
N_symb = N_symb-Lpreamble-Lnoisestart;


%%% Generation of the sequence of symbols (I)
M = 2;
N_bits_by_symb = log2(M);
I = pammod(N_symb,M);	% "pammod" is a function I wrote for the PAM modulation
I(1) = 1;   			% Because the first bit of the Access Address is always the opposite of the last one of the preamble.
I = [preamble;I];
LI = length(I);


%%% Gaussian pulse-shape filter (g) generation
N_Ts = 3;                   % Number of symbol periods for the filter support
time = -N_Ts/2*Ts : Ts/reso : N_Ts/2*Ts;    % Filter support
Tb = Ts/N_bits_by_symb;     % Here, we will always have Tb=Ts since M=2
K = sqrt(log(2));           % Parameter for the filter amplitude, set to a value giving a unitary amplitude
g = g_GMSK(time,K,Tb,BTb);  % Filter generation, using my function "g_GMSK"
Lg = length(g);


%%% Convolution (convo) between the symbol sequence and the Gaussian filter
Iup = upsample(I, reso);
convIupg = conv(Iup,g);
tmin = (reso*(LI-1)+Lg+1)/2 - (reso*LI)/2;
tmax = (reso*(LI-1)+Lg+1)/2 + (reso*LI)/2;
convo = convIupg(tmin:tmax);
Lconvo = length(convo);


%%% Integration (integ) of the convolution
integ = zeros(Lconvo-1,1);
incr = zeros(Lconvo-1,1);
integ(1) = (convo(1)+convo(2))/2 * (Ts/reso);
for i = 2:Lconvo-1
    incr(i) = (convo(i)+convo(i+1))/2 * (Ts/reso);
    integ(i) = integ(i-1) + incr(i);
end

% In practice the values in ak are tabled (they can be precomputed).
% They are used for the CFO correction.
integ2 = integ(1+reso:end) - integ(1:end-reso);
a = integ2 * pi*m/Tb;
ak = a(reso:reso+2*reso-1);


%%% Phase (phi) generation
phi = pi*m/Tb * integ;
Lphi = length(phi);


%%% Transmitted signal (s) generation
s = exp(1j*phi);
Ls = length(s);






%% Reception


OSF = 8;                % Oversampling factor (number of samples by symbol period at RX)


BERcorr = zeros(Nsnr,1);
BER = zeros(Nsnr,1);
df_hat = zeros(Nsnr,1);
ind_trig = zeros(Nsnr,1);
SFindi = zeros(Nsnr,1);
% Loop for each different value of the noise variance
% Since the non-idealities mitigations are sized for a SNR value of 13 or
% 14 dB, the loop should only take those values.
% for i = 1:Nsnr
SNR = 13;
for i = SNR+1
    
    %%% AWGN (n) generation
    N = N0(i)*OSF;
    n = awgn(Lnoisestart*OSF+Ls, 1, N);     % My function awgn generates the noise
    
    
    %%% Received signal (r) generation
    r = [zeros(Lnoisestart*OSF, 1) ; s] + n;
    Lr = length(r);


        %%% CFO        
        df = dfppm/1e06 * f_c;
        time = (1/(r_symb*OSF) : 1/(r_symb*OSF) : (Lr)/(r_symb*OSF))';
        rcfo = exp(1j*2*pi*df*time).*r;
        
       
        %%% Filtering of the received CFO'ed signal
        fcut = 0.6e06;
        Nper = 8;
        [h0] = LPF('Hamming', fcut, OSF, Nper);
        Lh0 = length(h0);
        rconvh00 = conv(rcfo,h0);
        rconvh0 = rconvh00(1+Nper/2*OSF:end-Nper/2*OSF);
        rcfof = rconvh0;


        %%% Preamble detection
        % Hermitian product
        OSR = OSF;
        Hprod = rcfof(1+OSR:end).*conj(rcfof(1:end-OSR));
        % IIR filter
        IIRb = 1;
        IIRa = zeros(OSR+1,1);
        IIRa(1) = 1;
        IIRa(OSR+1) = 0.875;
        IIRHprod = filter(IIRb,IIRa,Hprod);
        absIIRHprod = abs(IIRHprod);
        LabsIIRHprod = length(absIIRHprod);
        % Search for preamble
        if (SNR == 13)
            thresh = 37.5;
        elseif (SNR == 14)
            thresh = 32.5;
        end
        p = 1;
        booly = false;
        while(booly == false)
            if (absIIRHprod(p) > thresh)
                booly = true;
            elseif (p >= LabsIIRHprod)
                booly = true;
                fprintf('Error, after 1000 iterations, there is still no value that has reached the threshold!');
            else
                p = p+1;
            end
        end
        ind_trig(i) = p;    % This is the index that triggered the detection.
        
		
        %%% Symbol timing recovery
        Noffset = 5;        % We need to wait for 5 periods after the 
                            % preamble detection in order to target the 
                            % max of the 6th preamble period.
        offset = Noffset*OSF;
        maxxtemp = absIIRHprod(ind_trig(i)+offset);
        indtemp = ind_trig(i)+offset;
        kk = 1;
        while (kk < 1*OSF)
            if (absIIRHprod(ind_trig(i)+offset+kk) > maxxtemp)
                maxxtemp = absIIRHprod(ind_trig(i)+offset+kk);
                indtemp = ind_trig(i)+offset+kk;
            end
            kk = kk+1;
        end
        SFindi(i) = indtemp;

       
        %%% CFO correction        
        df_hat_temp1 = 0;
        Nperiod = 6;
        b = 1;
        counter = 1;
        symbr = zeros(Nperiod*OSF,1);
        symbrold = zeros(Nperiod*OSF,1);
        dfh = zeros(Nperiod*OSF,1);
        dfh1 = zeros(Nperiod*OSF,1);
        demodsymbrrold = zeros(Nperiod*OSF,1);
        dk = zeros(Nperiod*OSF,1);
        akind = zeros(Nperiod*OSF,1);
        hermi_unnorm = zeros(Nperiod*OSF, 1);
        hermi = zeros(Nperiod*OSF, 1);
        while(b <= Nperiod*OSF)
            
            % The first sample used for the CFO correction is the one 1
            % period after the preamble detection.
            symbr(b) = rcfof(b+ind_trig(i)+OSF-1);
            symbrold(b) = rcfof(b+ind_trig(i)+OSF-1-OSF);
            
            hermi_unnorm(b) = symbr(b)*conj(symbrold(b));
            hermi(b) = hermi_unnorm(b)/abs(hermi_unnorm(b));
            dk(b) = imag(hermi(b));

            % Simulates the index used to fetch the tabled values from the
            % LUTs.
            akind(b) = mod(b,2*OSF);
            if (akind(b) == 0)
                akind(b) = 2*OSF;
            end
            
            sinak = sin(ak(akind(b)));
            cosak = cos(ak(akind(b)));
            dfh1(b) = (dk(b)-sinak)/cosak /(2*pi*Tb);   % CFO estimate
            
            df_hat_temp1 = df_hat_temp1+dfh1(b);
            dfh(b) = df_hat_temp1/counter;      % Averaging of the estimate
            counter = counter + 1;
                                            
            b = b+1;
            
        end

        df_hat(i) = dfh(end);   % Final CFO estimate.

        
        rcorr = exp(-1j*2*pi*df_hat(i)*time((Lnoisestart+Lpreamble)*OSF+1 : end)).*rcfo((Lnoisestart+Lpreamble)*OSF+1 : end);
        % Filtering of the corrected signal
        rcorrconvh0 = conv(rcorr,h0);
        rcorrconvh0 = rcorrconvh0(1+Nper/2*OSF:end-Nper/2*OSF);
        rcorrf = rcorrconvh0;
        Lrcorrf = length(rcorrf);  
        

        
    %%% Filtering of the received signal (to simulate the case without CFO)
    rconvh0 = conv(r,h0);
    rconvh0 = rconvh0(1+Nper/2*OSF:end-Nper/2*OSF);
    rf = rconvh0;
    Lrf = length(rf);
        
        
    %%% Demodulation
    
    % Index for the decimation prior to the demodulation
    indidown = mod(SFindi(i)+2,OSF);
    if (indidown == 0)
        indidown = OSF;
    end
    
    rcorrfdown = downsample(rcorrf(indidown:end),OSF/2);
    Lrcorrfdown = length(rcorrfdown);
    if (mod(Lrcorrfdown,2)==0)
        demodcorr = imag(rcorrfdown(2:2:end).*conj(rcorrfdown(1:2:end-1)));
    else
        demodcorr = imag(rcorrfdown(2:2:end-1).*conj(rcorrfdown(1:2:end-2)));
    end
    Ldemodcorr = length(demodcorr);
    
    
    rfdown = downsample(rf(indidown:end),OSF/2);
    Lrfdown = length(rfdown);
    if (mod(Lrfdown,2)==0)
        demod = imag(rfdown(2:2:end).*conj(rfdown(1:2:end-1)));
    else
        demod = imag(rfdown(2:2:end-1).*conj(rfdown(1:2:end-2)));
    end
    demod = demod(1+Lpreamble+Lnoisestart:end);
    Ldemod = length(demod);
    
    
    %%% Decision
    Icorr_hat = sign(demodcorr);
    BERcorr(i) = sum(I(1+Lpreamble:Lpreamble+Ldemodcorr)~=Icorr_hat)/N_symb;
    I_hat = sign(demod);
    BER(i) = sum(I(1+Lpreamble:Lpreamble+Ldemodcorr)~=I_hat)/N_symb;



end



% BER vs Eb/N0 plot
figure;
semilogy(Eb_N0_dB,BERcorr,'xg','LineWidth',1.5,'MarkerSize',8);
hold on;
semilogy(Eb_N0_dB,BER,'xk','LineWidth',1.5,'MarkerSize',8);
axis([0 Nsnr 10/N_symb 1]);
grid;
legend('After CFO correction', 'Without CFO');
title( strcat('BER -- CFO = ', num2str(df)) );
xlabel('E_{b}/N_{0} [dB]');
ylabel('BER');
