uwb_channel_doc.m

Threshold selection for UWB TOA estimation based on kurtosis

% modified S-V channel model evaluation
%%Written by Sun Xu, Kim Chee Wee, B. Kannan & Francois Chin on 22/02/2005
function [h]=uwb_channel_doc(num_channels,cm_num)
%parameter:

%num_channels          number of channel impulse responses to generate
%h                               finial channel generated 
%clear;
no_output_files = 1; % non-zero: avoids writing output files of continuous-time responses

%num_channels = 10; % number of channel impulse responses to generate

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 params 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
% 
% 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, [1 num_channels], mean_excess_delay*[1 1]);
% grid on
% title('Excess delay (nS)')
% xlabel('Channel number')
% 
% figure(3); clf; 
% plot([1:num_channels], RMS_delay, [1 num_channels], mean_RMS_delay*[1 1]);
% grid on
% title('RMS delay (nS)')
% xlabel('Channel number')
% 
% figure(4); clf; 
% plot([1:num_channels], num_sig_paths, [1 num_channels], mean_sig_paths*[1 1]);
% 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, [1 num_channels], mean_sig_e_paths*[1 1]);
% 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 first 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