Add new file

GFSK_FixedPoint.m

+%% GFSK Fixed-point implementation
+
+%% Transmission
+
+%%% Transmission/modulation parameters
+
+% fixed-point precision
+prec_op_TX = 32;           % precision to be used for the operations performed on the considered signal
+prec_cst = 8;           % precision for the constants
+
+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 = 1e04;          % Number of symbols used for the simulation
+reso = 8;               % Resolution factor (number of samples by symbol period)
+dfppm = 20;
+
+
+%%% 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
+preamble = fi([1;-1;1;-1;1;-1;1;-1], 1, prec_cst, prec_cst-1);
+Lpreamble = length(preamble);
+N_symb = N_symb-Lpreamble;
+
+Lnoisestart = 92;
+N_symb = N_symb-Lnoisestart;
+
+%%% Generation of the sequence of symbols (I)
+bits = fi(round(rand(N_symb,1)), 1, prec_cst, prec_cst-1,'ProductMode','SpecifyPrecision','ProductFractionLength',16-1,'ProductWordLength',16,'SumMode', 'SpecifyPrecision','SumFractionLength',16-1,'SumWordLength',16);
+I0 = bits-0.5;
+I1 = I0*2;
+I1(1) = fi(1, 1, prec_cst, prec_cst-1,'ProductMode','SpecifyPrecision','ProductFractionLength',16-1,'ProductWordLength',16,'SumMode', 'SpecifyPrecision','SumFractionLength',16-1,'SumWordLength',16);
+I1(2) = fi(1, 1, prec_cst, prec_cst-1,'ProductMode','SpecifyPrecision','ProductFractionLength',16-1,'ProductWordLength',16,'SumMode', 'SpecifyPrecision','SumFractionLength',16-1,'SumWordLength',16);
+I = fi([preamble;I1], 1, prec_cst, prec_cst-1,'ProductMode','SpecifyPrecision','ProductFractionLength',16-1,'ProductWordLength',16,'SumMode', 'SpecifyPrecision','SumFractionLength',16-1,'SumWordLength',16);
+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;                    % Here, we will always have Tb=Ts
+K = sqrt(log(2));           % Parameter for the filter amplitude
+g = g_GMSK(time,K,Tb,BTb);  % Filter generation, using my function "g_GMSK"
+Lg = length(g);
+
+% g_GMSK can be precomputed and tabled since it is constant, so it does not
+% have to be calculated in fixed-point, only the final values are
+% converted.
+gfi = fi(g,1,16,16-1,'ProductMode','SpecifyPrecision','ProductFractionLength',16-1,'ProductWordLength',16,'SumMode', 'SpecifyPrecision','SumFractionLength',16-1,'SumWordLength',16);
+
+
+%%% Convolution (convo) between the symbol sequence and the Gaussian filter
+tmin = (reso*(LI-1)+Lg+1)/2 - (reso*LI)/2;
+tmax = (reso*(LI-1)+Lg+1)/2 + (reso*LI)/2;
+Iup = upsample(I, reso);
+convfiIupgfi = convfi(Iup,gfi',16,prec_op_TX);
+convofi = convfiIupgfi(tmin:tmax);
+Lconvo = length(convofi);
+
+
+%%% Integration (integfi) of the convolution (trapezium method)
+% First, we divide everything by 256 as it is the max value the integral
+% can output considering a packet of +-300 bits. This way, saturation is
+% avoided.
+convofi = bitsra(convofi, 8);
+integfi = fi(zeros(Lconvo-1,1),1,32,32-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_TX-1,'ProductWordLength',prec_op_TX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_TX-1,'SumWordLength',prec_op_TX);
+integfi(1) = bitsra(convofi(1)+convofi(2), log2(2*reso));  % Division is performed via arithmetic shift right
+for i = 2:Lconvo-1
+    integfi(i) = integfi(i-1) + bitsra(convofi(i)+convofi(i+1), log2(2*reso));
+end
+% To avoid saturation resulting from multiplication by pi/2, the signal
+% should be divided by 2.
+integfi = bitsra(integfi, 1);
+
+
+%%% Phase (phi) generation
+fipi_2 = fi(pi/2, 1, prec_cst, prec_cst-2,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_TX-1,'ProductWordLength',prec_op_TX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_TX-1,'SumWordLength',prec_op_TX);
+phifi = fipi_2 * integfi;
+Lphifi = length(phifi);
+
+
+%%% Transmitted signal (s) generation
+s = exp(1j*double(phifi)*2^9);
+Ls = length(s);
+
+
+
+
+
+%% Reception
+
+
+% Fixed-point precision
+prec_op_RX = 16;
+precRX = 16;
+
+
+% Table of values for the CFO correction:
+ak = [-0.894221291841193;...
+      -0.823171345055248;...
+      -0.624515667660133;...
+      -0.335010157734548;...
+                       0;...
+       0.335010157734548;...
+       0.624515667660133;...
+       0.823171345055248;...
+       0.894221291841193;...
+       0.823171345055248;...
+       0.624515667660133;...
+       0.335010157734548;...
+                       0;...
+      -0.335010157734548;...
+      -0.624515667660133;...
+      -0.823171345055248];
+sinak = sin(ak);
+cosak = cos(ak);
+sinakfi = fi(sinak,1,prec_cst,prec_cst-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+cosakfi = fi(cosak,1,prec_cst,prec_cst-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+pow_cosak2pifi_m1 = fi(1./(cosak*2*pi),1,prec_cst,prec_cst-1);
+
+
+BERcorrfi = zeros(Nsnr,1);
+BERfi = zeros(Nsnr,1);
+df_hatfi = zeros(Nsnr,1);
+ind_trig = zeros(Nsnr,1);
+SFindi = zeros(Nsnr,1);
+% Loop for each different value of the noise variance
+% for i = 1:Nsnr
+SNR = 13;
+for i = SNR+1
+
+    
+    %%% AWGN (n) generation
+    % noise can remain in floating-point, it doesn't matter as it will be
+    % converted as soon as it will be added to the signal.
+    N = N0(i)*reso;
+    n = awgn(Lnoisestart*reso+Ls, 1, N);
+    
+    
+    %%% Received signal (r) generation
+    r = [zeros(Lnoisestart*reso, 1) ; s] + n;
+    Lr = length(r);
+
+    
+        %%% CFO        
+        df = dfppm/1e06 * f_c;
+        time = (1/(r_symb*reso) : 1/(r_symb*reso) : (Lr)/(r_symb*reso))';
+        rcfo = exp(1j*2*pi*df*time).*r;
+        rcfofi0 = fi(rcfo,1,precRX,precRX-4,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        % /!\ abs(rcfof) can go up to 13 (result from 100000 iterations 
+        % with 300 bits and SNR=13dB) => a normalization by 16 should be
+        % performed.
+        rcfofi = fi(bitsra(rcfofi0,4),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        Lrcfo = length(rcfo);
+        
+       
+        %%% Filtering of the received CFO'ed signal
+        % This is done in digital, so from now on signals must be in
+        % fixed-point. But the filter can be precomputed, so only its final
+        % coefficient values have to be converted.
+        % Moreover, as the filter coefficients are <= 1, there is no need
+        % for normalization. But one is done nonetheless as rcfo needs one.
+        fcut = 0.6e06;
+        Nper = 12;
+        h0 = LPF('Hamming', fcut, reso, Nper);
+        h0fi = fi(h0,1,prec_cst,prec_cst-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        rconvfih0fi0 = convfi(rcfofi,h0fi',prec_op_RX,prec_op_RX);
+        rconvfih0fi = rconvfih0fi0(1+Nper/2*reso:end-Nper/2*reso);
+        rficfof = rconvfih0fi;
+        
+        
+
+
+        %%% Preamble detection
+        % Hermitian product
+        OSR = reso;
+        Hprodfi0 = rficfof(1+OSR:end).*conj(rficfof(1:end-OSR));
+        % IIR filter
+        IIRb = 1;
+        IIRa = zeros(OSR+1,1);
+        IIRa(1) = 1;
+        IIRa(OSR+1) = 0.875;
+        impu = eye(500,1);
+        impuresp = filter(IIRb, IIRa, impu);
+        % impuresp can be precomputed and tabled, so only the final values
+        % need to be converted into fixed-point.
+        impurespfi = fi(impuresp,1,prec_cst,prec_cst-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        % /!\ abs(IIRHprod) can go up to 350 (result from 100000 iterations 
+        % with 300 bits and SNR=13dB) => a normalization by 512 should be
+        % performed, normalization of Hprodfi. But it is already 16*16 = 
+        % 256 times smaller due to the previous normalization, so it needs 
+        % an additional normalization factor of only 2 to reach 512!
+        Hprodfi = fi(bitsra(Hprodfi0,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        IIRHprodfi0 = convfi(Hprodfi,impurespfi,precRX,prec_op_RX);
+        IIRHprodfi = IIRHprodfi0(1:end-(length(impu)-1));
+        absIIRHprodfi = abs(IIRHprodfi);
+        LabsIIRHprodfi = length(absIIRHprodfi);
+
+        
+        % Search for preamble
+        if (SNR==13)
+            thresh = 37.5;
+        elseif (SNR==14)
+            thresh = 32.5;
+        end
+        % Due to the normalization, the comparison must be done wrt a
+        % threshold that has undergone the same normalization.
+        threshfi = fi(thresh/2^(4+4+1),1,prec_cst,prec_cst-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_cst-1,'ProductWordLength',prec_cst,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_cst-1,'SumWordLength',prec_cst);
+        p = 1;
+        booly = false;
+        while(booly == false)
+            if (absIIRHprodfi(p) > threshfi)
+                booly = true;
+            elseif (p >= LabsIIRHprodfi)
+                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;
+        
+        
+        %%% Symbol timing recovery
+        Noffset = 5;
+        offset = Noffset*reso;
+        maxxtemp = absIIRHprodfi(ind_trig(i)+offset);
+        indtemp = ind_trig(i)+offset;
+        kk = 1;
+        while (kk < 1*reso)
+            if (absIIRHprodfi(ind_trig(i)+offset+kk) > maxxtemp)
+                maxxtemp = absIIRHprodfi(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;
+        afikind = zeros(Nperiod*reso,1);
+        symbrfi = fi(zeros(Nperiod*reso,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        symbrfiold = fi(zeros(Nperiod*reso,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);        
+        hermifi_unnorm = fi(zeros(Nperiod*reso,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        hermifi = fi(zeros(Nperiod*reso,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        demodsymbrfirfiold = fi(zeros(Nperiod*reso,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        dkfi = fi(zeros(Nperiod*reso,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        dfhfi = fi(zeros(Nperiod*reso,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        dfh1fi0 = fi(zeros(Nperiod*reso,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        dfh1fi = fi(zeros(Nperiod*reso,1),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        while(b <= Nperiod*reso)
+            
+            symbrfi(b) = rficfof(b+ind_trig(i)+reso-1);
+            symbrfiold(b) = rficfof(b+ind_trig(i)+reso-1-reso);
+            
+            hermifi_unnorm(b) = symbrfi(b)*conj(symbrfiold(b));
+            % The division for the normalization is somehow problematic,
+            % as I don't know how to handle this in Matlab, I "cheat" by
+            % using a double value for the divider.
+            hermifi(b) = hermifi_unnorm(b)/double(abs(hermifi_unnorm(b)));
+            demodsymbrfirfiold(b) = imag(hermifi(b));
+
+            dkfi0 = demodsymbrfirfiold(b);
+            % Must be divided by 2 to avoid saturation when subtracting it
+            % with sinakfik
+            dkfi(b) = bitsra(dkfi0,1);
+            
+            % Simulated the index to fetch the values from the LUTs.
+            afikind(b) = mod(b,2*reso);
+            if (afikind(b) == 0)
+                afikind(b) = 2*reso;
+            end
+            
+            sinakfik0 = sinakfi(afikind(b));
+            % Must be divided by 2 to avoid saturation when subtracting it
+            % from dkfi            
+            sinakfik = bitsra(sinakfik0,1);
+
+            dfh1fi0(b) = (dkfi(b)-sinakfik)*pow_cosak2pifi_m1(afikind(b));
+            
+            % Here, I accumulate Nperiod*reso (48) values, so I have to
+            % normalize by this value beforehand (48->64).
+            dfh1fi(b) = bitsra(dfh1fi0(b),6);
+            df_hat_temp1 = df_hat_temp1+dfh1fi(b);
+            dfhfi(b) = df_hat_temp1/counter;
+            counter = counter + 1;
+                                            
+            b = b+1;
+            
+        end
+        
+
+        df_hatfi(i) = dfhfi(end);
+        normfact = 2^7/Tb;
+
+        % Correction is done in RF, and the denormalization is performed at
+        % the same time using the factor normfact.
+        rcorr = exp(-1j*2*pi*(double(df_hatfi(i))*normfact)*time((Lnoisestart+Lpreamble)*reso+1 : end)).*rcfo((Lnoisestart+Lpreamble)*reso+1 : end);
+        rcorrfi0 = fi(rcorr, 1, precRX);
+        % Looks like convfi(rcorrfi,h0fi) can go up to 10, so it seems
+        % appropriate to normalize by 16.
+        rcorrfi = fi(bitsra(rcorrfi0, 4),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+        rcorrficonvfih0fi0 = convfi(rcorrfi,h0fi',precRX,prec_op_RX);
+        rcorrficonvfih0fi = rcorrficonvfih0fi0(1+Nper/2*reso:end-Nper/2*reso);
+        rcorrfif = rcorrficonvfih0fi;
+               
+        
+    %%% Filtering of the received signal (for the CFO-free case)
+    rfi0 = fi(r, 1, precRX);
+    % Looks like convfi(rfi,h0fi) can go up to 12, so it seems appropriate
+    % to normalize by 16.
+    rfi = fi(bitsra(rfi0, 4),1,precRX,precRX-1,'ProductMode','SpecifyPrecision','ProductFractionLength',prec_op_RX-1,'ProductWordLength',prec_op_RX,'SumMode', 'SpecifyPrecision','SumFractionLength',prec_op_RX-1,'SumWordLength',prec_op_RX);
+    rficonvfih0fi0 = convfi(rfi,h0fi',precRX,prec_op_RX);
+    rficonvfih0fi = rficonvfih0fi0(1+Nper/2*reso:end-Nper/2*reso);
+    rfif = rficonvfih0fi;        
+    
+
+    %%% Demodulation
+    
+    indidown = mod(SFindi(i)+2,reso);
+    if (indidown == 0)
+        indidown = reso;
+    end       
+    
+    rcorrfifdown = downsample(rcorrfif(indidown:end),reso/2);
+    Lrcorrfifdown = length(rcorrfifdown);
+    if (mod(Lrcorrfifdown,2)==0)
+        demodficorr = imag(rcorrfifdown(2:2:end).*conj(rcorrfifdown(1:2:end-1)));
+    else
+        demodficorr = imag(rcorrfifdown(2:2:end-1).*conj(rcorrfifdown(1:2:end-2)));
+    end    
+    Ldemodficorr = length(demodficorr);
+
+
+    rfifdown = downsample(rfif(indidown:end),reso/2);
+    Lrfifdown = length(rfifdown);
+    if (mod(Lrfifdown,2)==0)
+        demodfi = imag(rfifdown(2:2:end).*conj(rfifdown(1:2:end-1)));
+    else
+        demodfi = imag(rfifdown(2:2:end-1).*conj(rfifdown(1:2:end-2)));
+    end
+    demodfi = demodfi(1+Lpreamble+Lnoisestart:end);
+    Ldemodfi = length(demodfi);
+    
+
+    %%% Decision
+    Icorrfi_hat = sign(demodficorr);
+    BERcorrfi(i) = sum(sign(I(1+Lpreamble:Lpreamble+Ldemodficorr))~=Icorrfi_hat)/N_symb;
+    Ifi_hat = sign(demodfi);
+    BERfi(i) = sum(sign(I(1+Lpreamble:Lpreamble+Ldemodfi))~=Ifi_hat)/N_symb;
+    
+
+end
+
+
+% BER vs Eb/N0 plot
+figure;
+semilogy(Eb_N0_dB,BERcorrfi,'xg','LineWidth',1.5,'MarkerSize',8);
+hold on;
+semilogy(Eb_N0_dB,BERfi,'xk','LineWidth',1.5,'MarkerSize',8);
+axis([0 Nsnr 10/N_symb 1]);
+grid;
+title( strcat('BER -- df = ', num2str(df)) );
+xlabel('E_{b}/N_{0} [dB]');
+ylabel('BER');