%%%%%%%%%%%%%%%%%%%%%% Underwater Channel Simulation %%%%%%%%%%%%%%%%%%%%%%%
clc, clear all, warning off, close all,
tstart= tic;
cd(fileparts(which('channel_simulator')));
path(genpath([pwd,'/MatlabFunctions']), path);
%% read data from prm file:
file_name= 'example_channel';
[parameters, descriptions] = textread([file_name, '.prm'], '%f %s');
%% Deterministic channel parameters:
h0= parameters(1); % surface height (depth) [m]
ht0= parameters(2); % TX height [m]
hr0= parameters(3); % RX height [m]
d0= parameters(4); % channel distance [m]
k= parameters(5); % spreading factor
c= parameters(6); % speed of sound [m/s]
c2= parameters(7); % speed of sound in bottom [m/s] (>1500 for hard, < 1500 for soft)
cut= parameters(8); % do not consider arrivals whose relative strength is below this level
fmin= parameters(9); % minimum frequency [Hz]
B= parameters(10); % bandwidth [Hz]
fmax=fmin+B;
df= parameters(11); % frequency resolution [Hz]
f_vec=(fmin:df:fmax).'; Lf=length(f_vec);
fc=(fmax+fmin)/2;
f0=fmin;
f_vec2=(f0-B/2:df:f0+B+B/2); dif_f=f_vec2-fc;
dt= parameters(12); % time resolution [seconds]
T_SS= parameters(13); % coherence time of the small-scale variations [seconds]
t_vec=(0:dt:T_SS); Lt=length(t_vec);
%% Small-scale channel parameters:
sig2s= parameters(14); % variance of SS surface variations (surfampl^2/2?)
sig2b= parameters(15); % variance of SS bottom variations
B_delp= parameters(16); % 3-dB width of the p.s.d. of intra-path delays (assumed constant for all paths)
Sp= parameters(17); % number of intra-paths (assumed constant for all paths)
mu_p= parameters(18); % mean of intra-path amplitudes (assumed constant for all paths)
nu_p= parameters(19); % variance of intra-path amplitudes (assumed constant for all paths)
%% Large-scale channel parameters:
T_tot= parameters(20); % total duration of the simulated signal [seconds]
N_LS=round(T_tot/T_SS);
t_tot_vec=(0:dt:T_tot); Lt_tot=length(t_tot_vec);
h_bnd= parameters(21:22).'; % range of surface height (LS realizations \in h+h_band)
ht_bnd= parameters(23:24).'; % range of transmitter height
hr_bnd= parameters(25:26).'; % range of receiver height
d_bnd= parameters(27:28).'; % range of channel distance
sig_h= parameters(29); % standard deviation of LS variations of surface height
sig_ht= parameters(30); % standard deviation of LS variations of transmitter height
sig_hr= parameters(31); % standard deviation of LS variations of receiver height
sig_d= parameters(32); % standard deviation of LS variations of distance height
a_AR= parameters(33); % AR parameter for generating L-S variations (constant for variables h, ht, hr, d)
%% Large-Scale (L-S) & Small-Scale (S-S) Simulation Methods:
% 1)'sim_dir': SIMplified L-S model & DIRect S-S model
% 2)'sim_seq': SIMplified L-S model & Statistically EQuivalent S-S model
% 3)'bel_dir': BELlhop L-S model & DIRect S-S model
% 4)'bel_seq': BELlhop L-S model & Statistically EQuivalent S-S model
switch parameters(34)
case 1, method='sim_dir';
case 2, method='sim_seq';
case 3, method='bel_dir';
case 4, method='bel_seq';
end
method_LS= lower(method(1:3));
method_SS= lower(method(5:7));
%% Doppler parameters:
[Dopp_params] = textread([file_name, '.dop'], '%f');
Dopp_params= reshape(Dopp_params, Lt_tot, 10);
% drifting:
vtd_tot= Dopp_params(:,1).';
theta_td_tot= Dopp_params(:,2).';
vrd_tot= Dopp_params(:,3).';
theta_rd_tot= Dopp_params(:,4).';
% vehicular:
vtv_tot= Dopp_params(:,5).';
theta_tv_tot= Dopp_params(:,6).';
vrv_tot= Dopp_params(:,7).';
theta_rv_tot= Dopp_params(:,8).';
% surface:
Aw_tot= Dopp_params(:,9).';
fw_tot= Dopp_params(:,10).';
%% Large-scale loop:
H_LS= zeros(Lf, Lt*N_LS);
del_h=0; del_ht=0; del_hr=0; del_d=0;
h=h0; ht=ht0; hr=hr0; d= d0; % initial values
adopp0=zeros(1,50);
for LScount= 1:N_LS
rndvec= randn(1, 4);
del_h= a_AR * del_h + sqrt(1-a_AR^2)*sig_h*rndvec(1);
if (del_h > h_bnd(2))||(del_h < h_bnd(1)),
del_h= del_h- 2*sqrt(1-a_AR^2)*sig_h*rndvec(1);
end
htemp=h;
h= h0+del_h;
del_ht= a_AR * del_ht + sqrt(1-a_AR^2)*sig_ht*rndvec(2);
if (del_ht > ht_bnd(2))||(del_ht < ht_bnd(1)),
del_ht= del_ht- 2*sqrt(1-a_AR^2)*sig_ht*rndvec(2);
end
httemp=ht;
ht= ht0+del_ht;
del_hr= a_AR * del_hr + sqrt(1-a_AR^2)*sig_hr*rndvec(3);
if (del_hr > hr_bnd(2))||(del_hr < hr_bnd(1)),
del_hr= del_hr- 2*sqrt(1-a_AR^2)*sig_hr*rndvec(3);
end
hrtemp=hr;
hr= hr0+del_hr;
del_d= a_AR * del_d + sqrt(1-a_AR^2)*sig_d*rndvec(4);
if (del_d > d_bnd(2))||(del_d < d_bnd(1)),
del_d= del_d- 2*sqrt(1-a_AR^2)*sig_d*rndvec(4);
end
dtemp=d;
d= d0+del_d;
% Find Large-scale model parameters:
switch method_LS
case 'sim'
[lmean,taumean,Gamma,theta,ns,nb,hp]= mpgeometry(h,h-ht,h-hr,d,fc,k,cut,c,c2);
case 'bel'
[lmean,taumean,theta,ns,nb,hp,p0_ind]=runBellhop(h,h-ht,h-hr,d,fc,c,c2);
lmean=[lmean(p0_ind), lmean(1:p0_ind-1), lmean(p0_ind+1:end)];
taumean=[taumean(p0_ind), taumean(1:p0_ind-1), taumean(p0_ind+1:end)];
theta=[theta(p0_ind), theta(1:p0_ind-1), theta(p0_ind+1:end)];
ns=[ns(p0_ind), ns(1:p0_ind-1), ns(p0_ind+1:end)];
nb=[nb(p0_ind), nb(1:p0_ind-1), nb(p0_ind+1:end)];
hp=[hp(p0_ind), hp(1:p0_ind-1), hp(p0_ind+1:end)];
end
% ignore paths with delays longer than allowed by frequency resolution:
lmean= lmean(taumean<1/df);
theta= theta(taumean<1/df);
ns= ns(taumean<1/df);
nb= nb(taumean<1/df);
hp= hp(taumean<1/df);
taumean= taumean(taumean<1/df);
P=length(lmean); % number of paths
% Reference path transfer function:
H0= 1./sqrt(lmean(1)^k* (10.^(absorption(f_vec/1000)/10000)).^lmean(1) );
H= hp(1)*repmat( exp(-1j*2*pi*f_vec*taumean(1)) , 1, Lt);
%% Find small-scale model parameters:
sig_delp= sqrt(1/c^2*((2*sin(theta)).^2).*(ns*sig2s+nb*sig2b));
rho_p= exp(-((2*pi*f_vec).^2) * (sig_delp.^2/2));
rho_p_bb= exp(-(2*pi*(f_vec-fmin).^2) * (sig_delp.^2/2));
Bp= ((2*pi*f_vec*sig_delp).^2).*B_delp;
sig_p= sqrt(.5*(mu_p.^2.*Sp.*(1-rho_p.^2)+Sp.*nu_p.^2));
%% Find doppler rates:
% drifting:
vtd= vtd_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
theta_td= theta_td_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
vrd= vrd_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
theta_rd= theta_rd_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
% vehicular:
vtv= vtv_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
theta_tv= theta_tv_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
vrv= vrv_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
theta_rv= theta_rv_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
% surface:
Aw= Aw_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
fw= fw_tot(1+(LScount-1)*(Lt-1):1+LScount*(Lt-1));
vw= 2*pi*fw.*Aw;
% first path doppler:
vdrift= vtd.*cos(theta(1)-theta_td)-vrd.*cos(theta(1)+theta_rd);
adrift= vdrift/c;
vvhcl= 0;
avhcl= vvhcl/c;
vsurf= 0;
asurf= vsurf/c;
adopp= adrift + avhcl + asurf*ns(1);
eff_adopp = adopp0(1)+cumsum(adopp);
Dopp= exp(1j*2*pi*f_vec*(eff_adopp*dt));
adopp0(1)=eff_adopp(end);
H= H.*Dopp;
%% small-scale simulation:
for p=2:P;
switch method_SS
case 'dir', % DIRectly generate gamma:
gamma=zeros(Lf, Lt);
for counti=1:Sp;
gamma_pi= mu_p+nu_p*randn(1,Lt);
gamma_pi=repmat(gamma_pi, Lf, 1)*Sp;
deltau_pi=zeros(Lf, Lt);
w_delpi= sig_delp(p)*sqrt(1-exp(-1*pi*B_delp*dt)^2)*randn(1, 2*Lt);
temp_deltau_pi=filter(1, [1, -exp(-pi*B_delp*dt)], w_delpi);
for countf= 1:Lf
deltau_pi(countf, :)= temp_deltau_pi(Lt+1:end);
end
gamma=gamma+ gamma_pi.*exp(-1j*2*pi*repmat(f_vec, 1,Lt).*deltau_pi);
end
case 'seq', % use the Statistically EQuivalent model:
if sig_delp(p)*f0 < .4,
fprintf(['\n Warning! sig_delp*f0 < 0.4 for path ' num2str(p),' of L-S realization ', num2str(LScount),'. The statistically equivalent model i
评论0