% -------------
% Evaluation of Bit Synch. Loop Bandwidth
% Lorenzo Simone

clear all

% Simulation settings

CN0=44; % start C/NO (dBHz)
Bif=12.0e6; % IF bandwidth
F1=8.125e6; % reference frequency
adinp=7; % A/D input power in dBm
GAD=7.59; % A/D gain (quant/volt)
NMAX=200; % Num. of different CINO analyzed
fad=4*F1; % A/D clock frequency (Hz)
fad_2=fad/2;
fo=504000; % offset-frequency
fsub=131072; % subcarrier frequency (Pz)
BitRate=16384;
Pt=.5; % transition density
mod_idx=1.4; % sub-carrier modulation index (rad)
RxLoss=.5; % RxIloss due to short-loop implementation

N=fad/(2*fsub) % number of samples per subcarrier cycle 
Ns=1; % number of analysed subcarrier cycles

t=1/fad.2:1/fad_2:(Ns*N/fad-2); % Ns subcarrier cycles [r,cl=size(t);

% --- AGC WORD CHAIN A
AGCref=[545 540 535 528 520 507 499 486 466 447 426 401 376 346 325 296 271
250 225 205 185 168 155 140 130 120 115 110 105];
EbNoref=-.5:1:27.5;

for i_CN0=1:NMAX

        snr=CN0-10*log10(Bif); % Signal to Noise Ratio @ A/D input
        sigma2 = 10^((adinp - 10*log10(20))/10)/(1 + 10^(snr/10));
        sgm     = sqrt(sigma2);
        ampl = sqrt(2*(10^((adinp - 10*log10(20))/10) - sigma2));

        %---- SHORT LOOP SIMULATION ---
        Sl=ampl*cos(2*pi*fo*t+mod_idx*sin(2*pi*t*fsub)); 
        S2=ampl*sin(2*pi*fo*t+mod_idx*sin(2*pi*t*fsub));

        Rl=sign(cos(2*pi*t*fo)); 
        R2=sign(sin(2*pi*t*fo));

        Ik=S1.*R1+R2.*S2; 
        Qk=S2.*R1-S1.*R2;

        %---- SUB-CARRIER correlator output ----
        inph_lo=sign(sin(2*pi*t*fsub)); % in-phase subcar. ref.
        Is=GAD*Qk.*inph_lo; % subcarrier in-phase samples
        mean_Is=4*mean(Is);% subcarrier in-phase samples average @ correlator output
        A_bit=mean_Is;  % amplitude of the base-band samples

        % ------- AGC Look-Up Table ---------
        Sc_N0=CN0+10*log10(2*besselj(l,mod_idx)^2); 
        R=Sc_N0-10*log10(BitRate); 
        AGC_word(i_CN0) = interpl(EbNoref,AGCref,R);

        C_NO(i_CN0)=CN0; 
        A(i_CN0)=A_bit;

        CN0=CN0+0.1;

        end % (for i_CN0=1:NMAX)

% Loop Detector vs CINO
% ---------------------
for i=I:NMAX
        SC_N0(i)=C_N0(i)+10*log10(2*besselj(1,mod_idx)^2)-RxLoss;
        R(i)=(10^(Sc_N0(i)/10))/BitRate;        % Eb/NO
        kd(i)=0.5*(2*Pt)*(1/(2*pi))*A(i)*erf(sqrt(R(i)))*(F1/2)/BitRate; % phase detector gain
        if (AGC_word(i)>650) K(i)=2^-4; end 
        if (AGC_word(i)<=650 & AGC_word(i)>600) K(i)=2^-5; end 
        if (AGC_word(i)<=600 & AGC_word(i)>445) K(i)=2^-6; end 
        if (AGC_word(i)<=445 & AGC_word(i)>260) K(i)=2^-7; end 
        if (AGC_word(i)<=260) K(i)=2^-8; end

% Loop Gain
% ---------
Fck=(F1/2)/(2^3); 
Ts=1/BitRate; 
Tnco=1/BitRate; 
Knco=2*pi*(Fck/2^20)*Tnco;

% Loop BW vs Eb/NO
% ---------
for i=1:NMAX
        Bl(i)=(Knco*kd(i)*K(i))/(4*Ts);
end

% Loop Bw VS Eb/No
% ---------------
figure (1)
Eb_NO=Sc_N0-10.*log10(BitRate)+RxLoss;
plot(Eb_N0,Bl,'k');
x1abel ('E_b/N_0 (dB)')
ylabel ('BitSync. Loop Bandwidth (Hz)')
title (['Mod.Index = ',num2str(mod_idx), ' rad --- Bit Rate = ',num2str(BitRate),' bps'])
grid
xy=axis;
axis([min(Eb_N0) max(Eb_N0) xy(3) xy(4)])

% AGCword vs Eb/N0
% ---------------
figure (2)
plot(Eb_N0,AGC_word,'k');
xlabel ('Eb/N-0 (dB)')
ylabel ('AGC_w_o_r_d (quants)')
title (['Mod.Index = ',num2str(mod-idx),' rad --- Bit Rate = ',num2str(BitRate),' bps'])
grid
xy=axis;
axis([min(Eb_N0) max(Eb_N0) xy(3) xy(4)])

% AGCword vs C/NO
% ---------------
figure (3)
plot(C_N0,AGC_word,'k');
x1abel PCIN 0 (dB)I)
ylabel ('AGC_w_o_r_d (quants)')
title (['Mod.Index = ',num2str(mcd-idx),' rad --- Bit Rate = ',num2str(BitRate),' bps'])
grid
xy=axis;
axis(fmin(C_N0) max(C_N0) xy(3)-,xy(4)])

% Amplitude vs Eb/NO
% -----------------
figure (4) 
plot (Eb_N0,A,'k')
grid
title ('Mean @ Correlator Output') 
x1abel ('Eb/N_0 (dB)') 
ylabel ('Amplitude (Quants)')
xy=axis;
axis([min(Eb_N0) max(Eb_N0) xy(3) xy(4)])