sv_channel_model.m

Threshold selection for UWB TOA estimation based on kurtosis

% MATLAB program for "normal" UWB environments
% modifed S-V channel model evaluation
%
% Written by Sun Xu, Kim Chee Wee, B. Kannan & Francois Chin on 22/02/2005

clear all

no_output_files = 0; % non-zero: avoids writing output files of continuous-time responses
num_channels = 2; % number of channel impulse responses to generate: 100
randn('state',12); % initialize state of function for repeatability
rand('state',12); % initialize state of function for repeatability
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 );
fprintf(1,['Model Parameters\n', 'Lam = %.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,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,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);

%********************************************************************
% Testing and ploting
%********************************************************************
% ----------------------------channel energy-----------------------------------
channel_energy = sum(abs(h).^2);
h_len = length(h(:,1));
t = [0:(h_len-1)] * ts;     % for use in computing excess & RMS delays
excess_delay = zeros(1,num_channels);
RMS_delay = zeros(1,num_channels);
num_sig_paths = zeros(1,num_channels);
num_sig_e_paths = zeros(1,num_channels);

for k=1:num_channels  
    %----------------- determine excess delay and RMS delay-------------------
    sq_h = abs(h(:,k)).^2 / channel_energy(k);
    t_norm = t - t0(k); % remove the randomized arrival time of first cluster
    excess_delay(k) = t_norm * sq_h;
    RMS_delay(k) = sqrt( ((t_norm-excess_delay(k)).^2) * sq_h );
    
    % ------determine number of significant paths (paths within 10 dB from peak)----
    threshold_dB = -10; % dB
    temp_h = abs(h(:,k));
    temp_thresh = 10^(threshold_dB/20) * max(temp_h);
    num_sig_paths(k) = sum(temp_h > temp_thresh);
    
    % ------determine number of sig.paths (captures x % of energy in channel)-------
    x = 0.85;
    temp_sort = sort(temp_h.^2); % sorted in ascending order of energy
    cum_energy = cumsum(temp_sort(end:-1:1)); % cumulative energy
    index_e = min(find(cum_energy >= x * cum_energy(end)));
    num_sig_e_paths(k) = index_e;
end  % end of k=1:num_channels  

energy_mean = mean(10*log10(channel_energy));
energy_stddev = std(10*log10(channel_energy));
mean_excess_delay = mean(excess_delay);
mean_RMS_delay = mean(RMS_delay);
mean_sig_paths = mean(num_sig_paths);
mean_sig_e_paths = mean(num_sig_e_paths);

fprintf(1,'Model Characteristics\n');
fprintf(1,'Mean delays: excess (tau_m) = %.1f ns, RMS (tau_rms) = %1.f\n', mean_excess_delay, mean_RMS_delay);
fprintf(1,' # paths: NP_10dB = %.1f, NP_85%% = %.1f\n',mean_sig_paths, mean_sig_e_paths);
fprintf(1,'Channel energy: mean = %.1f dB, std deviation = %.1f dB\n',energy_mean, energy_stddev);
figure(1); clf; plot(t, abs(h)); grid on
title('Impulse response realizations')
xlabel('Time (nS)')

figure(2); clf; plot([1:num_channels],excess_delay,'b-',[1:num_channels],mean_excess_delay*ones(1,num_channels),'r-');
title('Excess delay (nS)');
grid on
xlabel('Channel number');

figure(3); clf; plot([1:num_channels], RMS_delay,'b-',[1:num_channels], mean_RMS_delay*ones(1,num_channels),'r-');
grid on
title('RMS delay (nS)')
xlabel('Channel number')

figure(4); clf; plot([1:num_channels], num_sig_paths,'b-', [1:num_channels], mean_sig_paths*ones(1,num_channels),'r-');
grid on
title('Number of significant paths within 10 dB of peak');
xlabel('Channel number')

figure(5); clf; plot([1:num_channels], num_sig_e_paths,'b-', [1:num_channels], mean_sig_e_paths*ones(1,num_channels),'r-');
grid on
title('Number of significant paths capturing > 85% energy')
xlabel('Channel number')

temp_average_power = sum((abs(h))'.*(abs(h))', 1)/num_channels;
temp_average_power = temp_average_power/max(temp_average_power);
average_decay_profile_dB = 10*log10(temp_average_power);
threshold_dB = -40;
above_threshold = find(average_decay_profile_dB > threshold_dB);
ave_t = t(above_threshold);
apdf_dB = average_decay_profile_dB(above_threshold);

figure(6); clf; plot(ave_t, apdf_dB); grid on
title('Average Power Decay Profile')
xlabel('Delay (nsec)')
ylabel('Average power (dB)')

if no_output_files,
return
end

%**************************************************************************
%Savinge the data
%**************************************************************************
%-------save continuous-time (time,value) pairs to files-------------------
save_fn = sprintf('cm%d_imr', cm_num);

% -----------A complete self-contained file for Matlab users----------------
save([save_fn '.mat'], 't', 'h','t_ct', 'h_ct', 't0', 'np', 'num_channels', 'cm_num');

% Three comma-delimited text files for non-Matlab users:
% File #1: cmX_imr_np.csv lists the number of paths in each realization
dlmwrite([save_fn 'np.csv'], np, ','); % number of paths

% File #2: cmX_imr_ct.csv can open with Excel
% n’th pair of columns contains the (time,value) pairs for the n’th realization save continous time data
th_ct = zeros(size(t_ct,1),3*size(t_ct,2));
th_ct(:,1:3:end) = t_ct; % time
th_ct(:,2:3:end) = abs(h_ct); % magnitude
th_ct(:,3:3:end) = angle(h_ct); % phase (radians)
fid = fopen([save_fn 'ct.csv'], 'w');

if fid < 0,
    error('unable to write .csv file for impulse response, file may be open in another application');
end

for k = 1:size(th_ct,1)
    fprintf(fid,'%.4f,%.6f,', th_ct(k,1:end-2));
    fprintf(fid,'%.4f,%.6f\r\n', th_ct(k,end-1:end)); %\r\n for Windoze end-of-line
end

fclose(fid);

% File #3: cmX_imr_dt.csv can open with Excel
% discrete channel impulse response magnitude and phase pair realization.
% the first column is time. phase is in radians

% -----------save discrete time data------------------
th = zeros(size(h,1),2*size(h,2)+1);
th(:,1) = t'; % the ?rst column is time scale
th(:,2:2:end) = abs(h); % even columns are magnitude
th(:,3:2:end) = angle(h); % odd columns are phase
fid = fopen([save_fn '_dt.csv'], 'w');

if fid < 0,
    error('unable to write .csv file for impulse response, file may be open in another application');
end

for k = 1:size(th,1)
    fprintf(fid,'%.4f,%.6f,', th(k,1:end-2));
    fprintf(fid,'%.4f,%.6f\r\n', th(k,end-1:end)); %\r\n for Windoze end-of-line
end

fclose(fid);
return; % end of program