clear;
% pd = makedist('uniform');
est_toa=zeros(1,100);
est_toa1=zeros(1,100);
est_toa2=zeros(1,100);
est_toa3=zeros(1,100);
est_toa5=zeros(1,100);
accumulator=0;
accumulator1=0;
accumulator2=0;
accumulator3=0;
accumulator5=0;
% rmse=zeros(1,7);
% rmse1=zeros(7,200);
% rmse2=zeros(7,200);
% rmse3=zeros(7,200);
rmse=zeros(1,7);
rmse1=zeros(1,7);
rmse2=zeros(1,7);
rmse3=zeros(7,200);
rmse5=zeros(1,7);
% rmse1=zeros(1,5);
% rmse2=zeros(1,5);
output=zeros(1,1600);
P=0;
% toa_firstpath=80e-9;
Pfa=0.05;
e=0.3;
f=0;
B=5e9%Hz
% Pfa1=0.1;
% nmax=0;
snr=[0 5 10 15 20 25 30]
%%------------------------transmit pulses
Tguard=80*10^-9; % guard time
Tf=2e-7; %time of every pulse
Tint=1e-9; % time of every integrator
o=0.4*10^-9; % pulse chip time
num_sample=Tf/Tint;
W=0
num_pulse=30
% phase=2*randi([0,1],1,num_pulse)-1;
% phase=0.25e-9*(2*randi([0,1],1,30)-1)
toa_firstpath=(80e-9);%+phase(1);
% time_hopping=randperm(200)-1;
% time_hopping=random('unid',200,1,30)-1;
t=0:0.125e-9:Tf*num_pulse;%3GHz
% transmit pulse sequence
for q=0:(num_pulse)-1
% w=((1-(4.*pi.*((t-(q*2e-7)-(8e-8)).^2)))./(o.^2)).*exp(-(2.*pi.*((t-(q*2e-7)-(8e-8)).^2))./(o.^2));
% w=(1-4*pi.*(((t-(q*2e-7)-(time_hopping(q+1)*Tint)-(8e-8)).^2)/(o.^2))).*exp(-2*pi*(((t-(q*2e-7)-(time_hopping(q+1)*Tint)-(8e-8)).^2)/(o.^2)))
w=(1-4*pi.*(((t-(q*2e-7)-(8e-8)).^2)/(o.^2))).*exp(-2*pi*(((t-(q*2e-7)-(8e-8)).^2)/(o.^2)))
%w=((1-(4.*pi.*((t-(q*2e-7)-(8e-8)-phase(q+1)).^2)))./(o.^2)).*exp(-(2.*pi.*((t-(q*2e-7)-(8e-8)-phase(q+1)).^2))./(o.^2));
% w=(phase(q+1))*(((1-(4.*pi.*((t-(q*2e-7)-8e-8).^2)))./(o.^2)).*exp(-(2.*pi.*((t-(q*2e-7)-8e-8).^2))./(o.^2)));
W=W+w
end
plot(t,W);
% stem(W)
%%%%%%%%%%%%%%%
% j=-1;
% for i=0:(num_pulse)-1
% j=j+1;
% for k=(1+1600*i):(1600+1600*i)
% if t(k)>=((80+200*j)*1e-9)&t(k)<=((80.05+200*j)*1e-9)
% W(k)=2e10;%6e18;
% else
% W(k)=0;
% end
% end
% end
% ---------------energy a pulse
fun = @(X) abs((1-(4.*pi.*(X.^2)/(o.^2))).*exp((-2.*pi.*(X.^2)/(o.^2)))).^2;
Ex=(0.5e7).*integral(fun,0,Tf);
%------------------------
for choice_snr=1:length(snr)
% for R=1:200
%-------------------------------------------
no_output_files = 1; % non-zero: avoids writing output files of continuous-time responses
num_channels = 100; % number of channel impulse responses to generate
randn('state'); % initialize state of function for repeatability
rand('state'); % initialize state of function for repeatability
cm_num = 4; % channel model number from 1 to 8
% get channel model params based on this channel model number
[Lam,lambda,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 );
fprintf(1,['Model Parameters\n' ...
' Lam = %.4f, lambda = %.4f, Lmean = %.4f, lambda_mode(FLAG) = %d\n' ...
' lambda_1 = %.4f, lambda_2 = %.4f, beta = %.4f\n' ...
' Gam = %.4f, gamma0 = %.4f, Kgamma = %.4f, sigma_cluster = %.4f\n' ...
' nlos(FLAG) = %d, gamma_rise = %.4f, gamma_1 = %.4f, chi = %.4f\n' ...
' m0 = %.4f, Km = %.4f, sigma_m0 = %.4f, sigma_Km = %.4f\n' ...
' sfading_mode(FLAG) = %d, m0_sp = %.4f, std_shdw = %.4f\n', ...
' kappa = %.4f, fc = %.4fGHz, fs = %.4fGHz\n'], ...
Lam,lambda,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);
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,lambda,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);
% 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
% correct for 1/N scaling imposed by decimation
% h = h * N; % normalized below..
% prepare to add the frequency dependency
K = 1; % K = Ko*Co^2/(4*pi)^2/d^n
% since the K is a constant, and the effect will be removed after
% normalization, so the K is set to be 1
h_len = length(h(:,1));
if (cm_num == 1|cm_num == 2| cm_num == 7|cm_num == 8|cm_num ==9)
[h]= uwb_sv_freq_depend_ct_15_4a(h,fc,fs,num_channels,kappa);
else
[h]= uwb_sv_freq_depend_ct_15_4a(h,fc,fs,num_channels,0);
end
%figure,stem(h)
for num_repeat=1:100
%received signal of multipath fading channel
received_signal=conv(h(:,num_repeat),W);
% t1=0:0.125e-9:(0.125e-9)*((length(received_signal)-1))
% figure,plot(t1,received_signal)
% figure;stem(received_signal)
% received_signal33=conv(h(:,num_repeat),W);
% add gaussian noise
noisy_receivedsignal = awgn(received_signal,snr(choice_snr),'measured');
% figure,plot(abs(noisy_receivedsignal))
% % noisy_receivedsignal=received_signal;
% % received_bandpass=noisy_receivedsignal;
% a7=xcorr(noisy_receivedsignal,W)
% %dd=length(W);
% aa7=find(abs(a7)>0.0001)
% for v=1:length(aa7)
% a7_new(v)=a7(aa7(v));
% end
% % a7_new=a7(1476:end)
% figure,plot(abs(a7_new))%(dd:end)))
bpFilt = designfilt('bandpassfir','FilterOrder',50, ...
'CutoffFrequency1',2e9,'CutoffFrequency2',7e9, ...
'SampleRate',15e9);
% pass received signal from bandpass filter
received_bandpass=filter(bpFilt,noisy_receivedsignal)
% received signal after bandpass filter
% out of energy detector
%%%%%%%%---------------------
% for i=1:1600
% for j=0:29
% P=P+received_bandpass(i+1600*j);
% end
% output(i)=P;
% P=0;
% end
% for j=1:200
% V(j)=(abs(output(1,1+(8*(j-1)))))^2+(abs(output(1,2+(8*(j-1)))))^2 ...
% +(abs(output(1,3+(8*(j-1)))))^2+(abs(output(1,4+(8*(j-1)))))^2 ...
% +(abs(output(1,5+(8*(j-1)))))^2+(abs(output(1,6+(8*(j-1)))))^2 ...
% +(abs(output(1,7+(8*(j-1)))))^2+(abs(output(1,8+(8*(j-1)))))^2;
% end
% S=sort(V);
%%%%%%%%---------------------
for i=1:num_pulse
if i==1
for j=1:200
% V(i,j)=(abs(received_bandpass(1,2+(3*(j-1)))))^2+(abs(received_bandpass(1,3+(3*(j-1)))))^2 ...
% +(abs(received_bandpass(1,4+(3*(j-1)))))^2;%+(abs(received_bandpass(1,5+(4*(j-1)))))^2;
V1(i,j)=(abs(received_bandpass(1,1+(8*(j-1)))))^2+(abs(received_bandpass(1,2+(8*(j-1)))))^2 ...
+(abs(received_bandpass(1,3+(8*(j-1)))))^2+(abs(received_bandpass(1,4+(8*(j-1)))))^2 ...
+(abs(received_bandpass(1,5+(8*(j-1)))))^2+(abs(received_bandpass(1,6+(8*(j-1)))))^2 ...
+(abs(received_bandpass(1,7+(8*(j-1)))))^2+(abs(received_bandpass(1,8+(8*(j-1)))))^2;
% +(abs(received_bandpass(1,9+(16*(j-1)))))^2+(abs(received_bandpass(1,10+(16*(j-1)))))^2 ...
% +(abs(received_bandpass(1,11+(16*(j-1)))))^2+(abs(received_bandpass(1,12+(16*(j-1)))))^2 ...
% +(abs(received_bandpass(1,13+(16*(j-1)))))^2+(abs(received_bandpass(1,14+(16*(j-1)))))^2 ...
% +(abs(received_bandpass(1,15+(1