diff --git a/GFSK_FixedPoint.m b/GFSK_FixedPoint.m
new file mode 100644
index 0000000000000000000000000000000000000000..bd2da8e4716bcb05216876c94847ebadb37b2255
--- /dev/null
+++ b/GFSK_FixedPoint.m
@@ -0,0 +1,380 @@
+%% 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');