%--------------channel modeling ------------------------------
clear all
num_channels =100; % number of channel impulse responses to generate: 100
cm_num =1; % channel model number from 1 to 8

% --------get channel model parameters based on this channel model number-----------
[Lam,Lmean,lambda_mode,lambda_1,lambda_2,beta,Gam,gamma_0,Kgamma, sigma_cluster,nlos,gamma_rise,gamma_1,chi,m0,Km,sigma_m0,sigma_Km, sfading_mode,m0_sp,std_shdw,kappa,fc,fs] = uwb_sv_params_15_4a( cm_num );
ts = 1/fs;     % sampling frequency

% --------get a bunch of realizations (impulse responses)---------------------
[h_ct,t_ct,t0,np] = uwb_sv_model_ct_15_4a(Lam,Lmean,lambda_mode,lambda_1, lambda_2,beta,Gam,gamma_0,Kgamma,sigma_cluster,nlos,gamma_rise,gamma_1,chi,m0,Km,sigma_m0,sigma_Km,sfading_mode,m0_sp,std_shdw,num_channels,ts);

% -------------change to complex baseband channel-----------------------------
h_ct_len = size(h_ct, 1);
phi = zeros(h_ct_len, 1);

for k = 1:num_channels
    phi = rand(h_ct_len, 1).*(2*pi);
    h_ct(:,k) = h_ct(:,k) .* exp(phi .* i);
end

% ----------now reduce continuous-time result to a discrete-time result---------
[hN,N] = uwb_sv_cnvrt_ct_15_4a( h_ct, t_ct, np, num_channels, ts );

if N > 1,
    h = resample(hN, 1, N); % decimate the columns of hN by factor N
else
    h = hN;
end

% ----------------add the frequency dependency---------------------------------
[h]= uwb_sv_freq_depend_ct_15_4a(h,fc,fs,num_channels,kappa);

%---------------------------inital-------------------
fc = 8e9;      % sampling frequency
Ts = 1e-9;
% B = 4e9;
%{
Tf=200e-9;
Tb=4e-9;
Intagtime=Ts*fc ;
Nb=1.5*Tf/Tb;
%}
Intagtime = Ts*fc;
Nktosis = 300;

%_-----------transmitter------------------
[bits,THcode,Stx,ref]=cp0201_transmitter_2PPM_TH;  %generate the signal;  
Delay = randint(1,1,200); % generate the Delay
DelayStx=[Delay ref];  % Add the Delay;

ebn0 = 4:2:32;
 for i = 1:length(ebn0)
    for count=1: num_channels
        Channel=abs(h(:,count)); 
        DelayStxChannel=conv(DelayStx,Channel');
%         numbits=1;
%         Eb = (1/numbits)*sum(DelayStxChannel.^2);     % measured energy per bit
%         EbNo = 10.^(ebn0(i)./10);          % Eb/No in linear unit     
%         No = Eb ./ EbNo;                % Unilateral spectral
%                                  % density
%         % noise=zeros(length(EbNo),length(DelayStxChannel));
%         % r =  zeros(length(EbNo),length(DelayStxChannel));                       
%         nstdv = sqrt(No./2);       % Standard deviation for the noise       
%         noise = nstdv .* randn(1,length(DelayStxChannel));
        %-----------awgn---------------
       spow = sum(DelayStxChannel.^2)/length(DelayStxChannel)*512;
       attn = 0.5*spow*10.^(-ebn0(i)/10);
       attn = randn(1,length(DelayStxChannel))*sqrt(attn);
       %---------------------rx-----------------
       r = DelayStxChannel + attn;    
       %compute the Kuritos
       Karray=r(1:Nktosis*Intagtime);
       Karray= Karray.^2;  % Integral 
       Karray=reshape(Karray,Intagtime,Nktosis);
       KarrayInt= abs(sum(Karray./fc));     
       KNktosis(count) = kurtosis(KarrayInt)-3;
       %  KNktosis(count)=kurtosis(KarrayInt)-3;
    end     
     K_KNktosis(i) = mean(KNktosis);
  end
 plot(ebn0,K_KNktosis,'g:o');
 grid on;
 title('Figure 1: The expected value of the kurtosis vs Eb/N0 for different channel models and block sizes in energy detection' );
 xlabel('EbN0 (dB)');
 ylabel('E(K)');