%=============================================================================%
%=                         TOA UWB ALGORITHM                                 =%
%=                                                                           =%
%= Programmed by Jasurbek Khodjaev                                           =%
%= Yeungnam University Mobile Communication Lab.                             =%
%= MCL 2006                                                                  =%
%=============================================================================%

clear all;
clc;

%--------------- Initialization -----------------------------------------
%-------------------------------------------

% Speed of Light
light_speed = 3e8;

% Coordinates of APs
AP = [0 0; 0 10; 10 0]; % in meters

% Number of Access Points (AP)
num_ap = length(AP);

% Tag's initial coordinate
% Tag = [5 4];



    
    
% Pulse shape
pulse_order = 1; % 0-Gaussian pulse, 1-First derivative of Gaussian pulse, 2 - Second derivative;

% Number of bits
num_bits = 1000;

% Pulse repetition interval, PRI
pri = 200e-9;

% The SNR range (in dB)
EbNo = -15;

fs = 10e9; %sample rate-10 times the highest frequency in GHz
ts = 1/fs; %sample period
t = [(-0.5E-9):ts:(0.5E-9)]; %vector with sample instants
t1 = .5E-9; %pulse width(0.5 nanoseconds)
%-------------------------------------------------------------------------


%----------------- PRI -------------------------------------------

A =-1;%positive value gives negative going monopulse;neg value gives
   %positive going monopulse
[y] = monocycle(fs, ts, t, t1, A, pulse_order); % Generate Gaussian pulse

n_pulse_pri = round(pri/ts);               % Sampling of PRI
sig = zeros(1,n_pulse_pri);    
sig(1:length(y)) = y;                 % One pulse in one PRI
%-----------------------------------------------------------------
shift_const = 40;

kkk = 0;
for EbNo = -15:5:5
    kkk = kkk + 1;
    
        % Tag's initial coordinate
    a = round(rand * 10);
    b = round(rand * 10);
    Tag = [a b];
    
% Distance calculation between each AP and the Tag, IDEAL case
for ii = 1:num_ap
    dist_ap_tag(ii) = dist_t(AP(ii,:), Tag);
    % Time from each AP to Tag
    time_ap_tag(ii) = dist_ap_tag(ii)/light_speed;
end




    %++++++++++++++++ TRANSMISSION +++++++++++++++++++++++++++++++
            noise_var   = 10^(-EbNo/10);

%     y=y/std([y zeros(1,4000-length(y))]);
    for jj = 1:num_bits
            % From TAG to AP1
            del_sample_ap_tag1 = round(time_ap_tag(1)/ts);
            xx1 = zeros(1,del_sample_ap_tag1);
            del_sig_ap1_tag(jj,:) = [ sig(1:end-length(xx1)) xx1];
            h1 = uwb_channel(dist_ap_tag(1));
            h1=[zeros(1,del_sample_ap_tag1) h1];
            h1=h1+randn(1,length(h1)) .* sqrt(noise_var)*std(h1);
            
            conv_data1 = conv(del_sig_ap1_tag(jj,:), h1);
            ap1_tag_chan(jj, :) = conv_data1(1:length(sig));

            % From TAG to AP2
            del_sample_ap_tag2 = round(time_ap_tag(2)/ts);
            xx2 = zeros(1,del_sample_ap_tag2);
            del_sig_ap2_tag(jj,:) = [ sig(1:end-length(xx2)) xx2];
            h2 = uwb_channel(dist_ap_tag(2));
            h2=[zeros(1,del_sample_ap_tag2) h2];
            h2=h2+randn(1,length(h2)) .* sqrt(noise_var)*std(h2);
            
            conv_data2 = conv(del_sig_ap2_tag(jj,:), h2);
            ap2_tag_chan(jj, :) = conv_data2(1:length(sig));
            
            % From TAG to AP3
            del_sample_ap_tag3 = round(time_ap_tag(3)/ts);
            xx3 = zeros(1,del_sample_ap_tag3);
            del_sig_ap3_tag(jj,:) = [ sig(1:end-length(xx3)) xx3];
            h3 = uwb_channel(dist_ap_tag(3));
            h3=[zeros(1,del_sample_ap_tag3) h3];
            h3=h3+randn(1,length(h3)) .* sqrt(noise_var)*std(h3);
            
            conv_data3 = conv(del_sig_ap3_tag(jj,:), h3);
            ap3_tag_chan(jj, :) = conv_data3(1:length(sig));
           [kkk jj]
%     end

    %-------------------------------------------------------
    % Additive White Gaussian Noise (AWGN) Channel ---------
%     noise_var   = 10^(-EbNo/10);
%     for jj = 1:num_bits
        ap1_tag_chan_wgn(jj,:) = ap1_tag_chan(jj,:) + 0*std([y zeros(1,2000-length(y))])* randn(1,length(ap1_tag_chan(jj,:))) .* sqrt(noise_var);
        ap2_tag_chan_wgn(jj,:) = ap2_tag_chan(jj,:) + 0*std([y zeros(1,2000-length(y))])*randn(1,length(ap2_tag_chan(jj,:))) .* sqrt(noise_var);
        ap3_tag_chan_wgn(jj,:) = ap3_tag_chan(jj,:) + 0*std([y zeros(1,2000-length(y))])*randn(1,length(ap3_tag_chan(jj,:))) .* sqrt(noise_var);
        
%         ap1_tag_chan_wgn(jj,:) = ap1_tag_chan(jj,:)/std(ap1_tag_chan(jj,:)) + randn(1,length(ap1_tag_chan(jj,:))) .* sqrt(noise_var);
%         ap2_tag_chan_wgn(jj,:) = ap2_tag_chan(jj,:)/std(ap2_tag_chan(jj,:)) + randn(1,length(ap2_tag_chan(jj,:))) .* sqrt(noise_var);
%         ap3_tag_chan_wgn(jj,:) = ap3_tag_chan(jj,:)/std(ap3_tag_chan(jj,:)) + randn(1,length(ap3_tag_chan(jj,:))) .* sqrt(noise_var);
%         ap1_tag_chan_wgn(jj,:) = cp0801_Gnoise1(ap1_tag_chan(jj,:), EbNo, 1);
%         ap2_tag_chan_wgn(jj,:) = cp0801_Gnoise1(ap2_tag_chan(jj,:), EbNo, 1);
%         ap3_tag_chan_wgn(jj,:) = cp0801_Gnoise1(ap3_tag_chan(jj,:), EbNo, 1);
%         
    end
    
% -------------------------------------------------------

%     received_signl_ap1 = sum(ap1_tag_chan_wgn)/num_bits;
%     xc = cp0804_corrsyn(received_signl_ap1, sig, fs);
% %     cp0804_corrsyn(signal,template,fc)
%     [a,delay1]=max(xc);
%     TOA_1 = (length(sig) - delay1) * ts;
%     
%     %-----------------------------------
%     
%     received_signl_ap2 = sum(ap2_tag_chan_wgn)/num_bits;
%     xc = cp0804_corrsyn(received_signl_ap2, sig, fs);
%     [a,delay2]=max(xc);
%     TOA_2 = (length(sig) - delay2) * ts;
%     
%     %-----------------------------------
%     
%     received_signl_ap3 = sum(ap2_tag_chan_wgn)/num_bits;
%     xc = cp0804_corrsyn(received_signl_ap3, sig, fs);
%     [a,delay3]=max(xc);
%     TOA_3 = (length(sig) - delay3) * ts;
    
    

    %------------------- AP1 receiver ------------------------------
    % Correlator
    received_signl_ap1 = sum(ap1_tag_chan_wgn)/num_bits;
    xc1 = xcorr(y, received_signl_ap1);
    [a,delay1]=max(xc1(1:5:2000));
    TOA_1 = (length(sig)/5 - delay1) * ts*5;
%     TOA_1 = delay1 * ts;

    %------------------- AP2 receiver ------------------------------
    % Correlator
    received_signl_ap2 = sum(ap2_tag_chan_wgn)/num_bits;
    xc2 = xcorr(y, received_signl_ap2);
    [a,delay2]=max(xc2(1:5:2000));
    TOA_2 = (length(sig)/5 - delay2) * ts*5;
%     TOA_2 =  delay2 * ts;

    %------------------- AP3 receiver ------------------------------
    % Correlator
    received_signl_ap3 = sum(ap3_tag_chan_wgn)/num_bits;
    xc3 = xcorr(y, received_signl_ap3);
    [a,delay3] = max(xc3(1:5:2000));
    TOA_3 = (length(sig)/5 - delay3) * ts*5;
    
%     TOA_3 =  delay3 * ts;
    %---------------------------------------------------------------
%     time_ap_tag = time_ap_tag;
    time_dur = [TOA_1 TOA_2 TOA_3];
    [ex,ey]=TOA_LS(0,AP(:,1),AP(:,2),time_dur*3e8)
    toa_error(1,kkk) = toa(AP, Tag, time_dur, light_speed);
    kkk
end
% 

    EbNo = -15:5:5;
    figure,plot(EbNo , toa_error);