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