Add new file

GFSK_FloatingPoint

+%% 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);
+    
+    
+    %%% 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');