ber_sim_awgn_multipath_data_154a.m

一个UWB仿真程序包

%
%
function [num_err,ber_sim,timing,chunk] = ...
    ber_sim_awgn_multipath_data_154a(FS_CONT,DATA_LENGTH,RSCODE, ...
                                     Nu, LAMBDA, Npacket, ... 
                                     TC,BURST_LENGTH,SYMBOL_LENGTH, ...
                                     CH_MODEL,CH_ATT_THLD_DB, channel_est, ...
                                     snr_dB,B, ...
                                     T_INT,N_INT, ...
                                     timing_offset, ...
                                     RESULT_FILE)

%
% Necessary path
addpath ../../matlab/util ../../matlab/rng ../../matlab/signal ...
    ../../matlab/receiver ../../matlab/data_signal ...
    ../../matlab/reed_solomon ...
    ../../src/channel ../../src/scrambling_timehopping ...
    ../../matlab/preamble_signal ../../matlab/mui;


addpath ../../../UWB/sim/matlab/toolbox_comm;
addpath ../../../UWB/sim/matlab/toolbox_signal;
addpath ../../../UWB/sim/matlab/toolbox_stats;


fprintf('FS_CONT = %2.0e, DATA_LENGTH = %d, Npacket = %d\n', ...
        FS_CONT,DATA_LENGTH,Npacket);
fprintf('Number of users = %d, Rate = %.2f [pkts/sec]\n',Nu,LAMBDA);
fprintf('TC = %2.1e, BURST_LENGTH = %d, SYMBOL_LENGTH = %d\n', ...
        TC,BURST_LENGTH,SYMBOL_LENGTH);
fprintf('CH_MODEL = %d, CH_ATT_THLD_DB = %d\n',CH_MODEL,CH_ATT_THLD_DB);
if channel_est == 1
  fprintf('Do channel estimation\n');
else
  fprintf('No channel estimation\n');
  fprintf('Length of the mask = %e [s]\n',channel_est*T_INT*N_INT);
end
fprintf('snr  = %d [dB], B = %2.1e\n',snr_dB,B);
fprintf('T_INT = %2.1e, N_INT = %d\n',T_INT,N_INT);
fprintf('timing_offset = %2.1e\n',timing_offset);

if RSCODE == 1
  fullblock = floor(DATA_LENGTH/330);
  PACKET_LENGTH = fullblock*378 + (DATA_LENGTH-fullblock*330) + 48;
  fprintf('RS coding enabled, ');
else
  PACKET_LENGTH = DATA_LENGTH;
  fprintf('No RS coding, ');
end
fprintf('PACKET_LENGTH = %d\n',PACKET_LENGTH);


% Random number generator initialization
rand('seed',prod(clock));
randn('seed',prod(clock));


% $$$ % Simulation parameters
% $$$ FS_CONT     = 12e9;
% $$$ PACKET_LENGTH = 127*8;
% $$$ Npacket      = 16;
% $$$ % Physical layer
% $$$ TC            = 2e-9; % Chip length in seconds
% $$$ BURST_LENGTH  = 4; % LPRF, in chip
% $$$ %BURST_LENGTH  = 16; % HPRF, in chip
% $$$ SYMBOL_LENGTH = 512; % In chip
% $$$ % Multipath channel
% $$$ CH_MODEL       = 1;
% $$$ CH_ATT_THLD_DB = 30; %rays with bigger attenuation are considered 0
% $$$ CH_RNG_SEED    = sum(100*clock);
% $$$ % Bandwidth of the noise
% $$$ snr_dB = 10;
% $$$ snr    = 10.^(snr_dB/10);
% $$$ Eb     = 1;
% $$$ s2     = Eb ./ (2*snr);
% $$$ B      = 2e9;
% $$$ % Receiver
% $$$ T_INT = 1e-9;
% $$$ N_INT = 4;
% $$$ % Result from the timing synchronization and from the channel
% $$$ % estimation
% $$$ channel_mask       = [];
% $$$ % $$$ channel_mask       = ones(1,16);
% $$$ % $$$ channel_mask_unit  = 8e-9; % Duration in seconds of a channel mask
% $$$ % $$$                            % sample
% $$$ channel_mask_unit  = T_INT*N_INT;
% $$$ timing_offset      = 0;  % In seconds
% $$$ timing_offset_samples = timing_offset*FS_CONT;

% Variance of the noise
snr    = 10.^(snr_dB/10);
Eb     = 1;
s2     = Eb ./ (2*snr);
% Seed for the channel impulse response generator
CH_RNG_SEED    = sum(100*clock);

%group all the parameters corresponding to the channel
channel_p = group_channel_params(CH_RNG_SEED,CH_MODEL,CH_ATT_THLD_DB,0.5,1e-9);


% Channel mask if any
% channel_mask       = [];
% $$$ channel_mask       = ones(1,16);
% $$$ channel_mask_unit  = 8e-9; % Duration in seconds of a channel mask
% $$$                            % sample
% $$$ channel_mask_unit  = T_INT*N_INT;
% Timing-error stemming from the timing acquisition
timing_offset_samples = timing_offset*FS_CONT;

t0 = clock;

% Pre-construction of the data signal of interest
samples = get_samples_values(TC,FS_CONT, ...
                                BURST_LENGTH, ...
                                SYMBOL_LENGTH, ...
                                PACKET_LENGTH);

if Nu > 1
  % Draw the time of arrival of the packets
  PKTS_PER_RUN = Npacket;
  WARMUP = 500;
  L = 64; %spreading factor
  G = 4;
  CODE_ID = 5;
  % Get preamble values
  [uoi_pre sfd_code NSYNC NSFD] = generate_preamble(L,CODE_ID);
  code_len = length(len31_preamble_code(CODE_ID));

  ack_gap = ((NSYNC + NSFD)*code_len*L*TC) + (0*SYMBOL_LENGTH*TC);
  % Get the arrival time of all the packets
  [timing] = get_packet_timing(Nu,PKTS_PER_RUN,LAMBDA,ack_gap, ...
                                  WARMUP,NSYNC,NSFD,TC,L,G, ...
                                  PACKET_LENGTH,SYMBOL_LENGTH, ...
                                  code_len,0);

  % Best we can do is a precision of 1/FS_CONT -> convert cont. times to this
  % precision
  timing.packet_len   = round(timing.packet_len*FS_CONT);
  timing.preamble_len = round(timing.preamble_len*FS_CONT);
  timing.payload_len  = round(timing.payload_len*FS_CONT);
  timing.packets(:,1) = round(timing.packets(:,1)*FS_CONT);

  % Add some random samples at the beginning so that we are not always well
  % aligned with a chip -> effect of arbitrary sampling start
  %timing.startTimes = timing.startTimes - ...
  %    floor(samples.chip*rand(1,length(timing.startTimes)));

  % Create the simulation chunks
  chunk.idx    = find(timing.packets(:,2) == 1);
  chunk.length = length(chunk.idx);
  chunk.start  = timing.packets(chunk.idx,1)+timing.preamble_len;
  chunk.stop   = timing.packets(chunk.idx,1)+timing.packet_len;
% $$$   for l = 1:chunk.length
% $$$     fprintf('Chunk (%d): starts at %12.2f, ends %12.2f\n', ...
% $$$             l,chunk.start(l),chunk.stop(l));
% $$$   end

  chunk.mui_packet = zeros(1,chunk.length);
  for l = 1:chunk.length
    chunk.mui_delay(l).delay = [];
  end
  % Iterate through all the packets that are not from the UOI
  mui_idx = find(timing.packets(:,2) ~= 1);
  for k = 1:length(mui_idx)
    fprintf('MUI packet %d (%d) starts at %12.2f\n', ...
            mui_idx(k),timing.packets(mui_idx(k),2),timing.packets(k,1));
    for l = 1:chunk.length
% $$$       if (timing.packets(mui_idx(k),1) > chunk.start(l) ...
% $$$           & timing.packets(mui_idx(k),1) < chunk.stop(l)) ...
% $$$             | (timing.packets(mui_idx(k),1)+timing.packet_len > chunk.start(l) ...
% $$$                & timing.packets(mui_idx(k),1)+timing.packet_len < chunk.stop(l))
      mui_delay = chunk.start(l)-timing.packets(mui_idx(k),1);
      if (sign(mui_delay) == 1 & mui_delay<timing.packet_len) ...
            | (sign(mui_delay) == -1 & -mui_delay<timing.payload_len)
        % This MUI packet overlaps possibly
        fprintf('MUI packet %d overlaps with chunk %d, ',mui_idx(k),l);
        % This is a debugging test
        if sign(mui_delay) == 1
          if (timing.packets(mui_idx(k),1)+timing.packet_len > ...
                chunk.stop(l)+timing.preamble_len)
            fprintf('Alert!!!!\n');
	    exit;
          end
        elseif sign(mui_delay) == -1
          if timing.packets(mui_idx(k),1) < chunk.start(l)
            fprintf('Alert!!!!\n');
	    exit;
          end
        end
        chunk.mui_packet(l) = chunk.mui_packet(l) + 1;
        chunk.mui_delay(l).delay = ...
            [chunk.mui_delay(l).delay mui_delay];
        fprintf('delay %12.2f\n',chunk.mui_delay(l).delay(end));
      end
    end
  end
  % Print the number of overlapping MUI packets
else
  timing = [];
  chunk  = [];
end

[ep] = get_ep(PACKET_LENGTH, ...
              SYMBOL_LENGTH, ...
              BURST_LENGTH, ...
              samples);

symbol_start = get_symbol_start(PACKET_LENGTH,samples);


[burst, compound_channel, burstlength] ...
    = init_data_signal(FS_CONT, PACKET_LENGTH, TC,BURST_LENGTH,SYMBOL_LENGTH, ...
                       ep.ths, ep.scrambling, samples, ...
                       channel_p);

% Parameters for the generation of the Gaussian noise
ALPHA = PACKET_LENGTH/2;
Ts    = TC * SYMBOL_LENGTH;
Nt    = Ts * FS_CONT;

% Get the data bit (coded or uncoded)
[tx_bit,data_bit] =  get_tx_bit(DATA_LENGTH,RSCODE,ep.interleaver_seed);

% Some stat
overlap.do       = zeros(1,Npacket);
overlap.preamble = zeros(1,Npacket);
overlap.payload  = zeros(1,Npacket);
overlap.err      = zeros(1,Npacket);
overlap.ber      = zeros(1,Npacket);

for l = 1:Npacket
  fprintf('\n');
  fprintf('> Packet %d, ',l);

  if channel_est == 1
    % Fake estimation of the channel mask
    L             = 64;
    CODE_ID       = 5;
    preamble_code = abs(len31_preamble_code(CODE_ID));
    channel_mask  = fake_estimate_channel(TC,FS_CONT,samples,...
					  L,preamble_code,CODE_ID,...
					  compound_channel,...
					  T_INT,N_INT,...
					  s2,B);
    channel_mask  = channel_mask';
    channel_mask_unit  = T_INT*N_INT;
  else
    channel_mask = ones(1,channel_est);
    channel_mask_unit  = T_INT*N_INT;
  end
  channel_mask_length(l) = length(find(channel_mask)>0);

  % Get the data signal of interest
  [rx_data, signal_start] ... 
      = get_data_signal(PACKET_LENGTH, ...
                        samples.burst, samples.symbol, burstlength, ...
                        burst,ep.ths,symbol_start,tx_bit);
  %size(rx_data)
  if Nu > 1

    % Add the signal generated by the MUI packet to the the received signal
    for k = 1:chunk.mui_packet(l)
      fprintf('Chunk %d: generate Mui packet %d\n',l,k);

      % Stat
      overlap.do(l) = overlap.do(l) + 1;

% $$$     packet = create_packet(start_samples,stop_samples,compound_channel,...
% $$$     ch_len,chip_span,tx_bit,ep.ths,ep.scrambling,symbol_start,CODE_ID,L,TC,...
% $$$     FS_CONT,DATA_LENGTH,RSCODE,PACKET_LENGTH,BURST_LENGTH,SYMBOL_LENGTH,...
% $$$     samples,PLOT_DEBUG)
      if (chunk.mui_delay(l).delay(k) > 0)
        fprintf('MUI packet starts before the data packet\n');
        fprintf('Delay = %12.2f (%.2f)\n',chunk.mui_delay(l).delay(k),...
                1-(chunk.mui_delay(l).delay(k)/timing.payload_len));

        start_samples = chunk.mui_delay(l).delay(k);
        stop_samples  = min(timing.packet_len,start_samples+timing.payload_len-1);

	% Stat
	% Do we overlap with the preamble, or the data part. or both?
        % Valid only for Nu == 2
        if timing.packet_len > start_samples + timing.payload_len
          fprintf('[%d] Full overlap\n',overlap.do(l));
        end
	if start_samples <= timing.preamble_len
	  fprintf('[%d] Overlap with preamble and data\n',overlap.do(l));

	  ov = (timing.preamble_len-start_samples+1)/timing.payload_len;
	  overlap.preamble(l) = overlap.preamble(l)+ov;

	  ov = (stop_samples - timing.preamble_len+1)/timing.payload_len;
	  overlap.payload(l) = overlap.payload(l)+ov;
	else
	  fprintf('[%d] Overlap with data only\n',overlap.do(l));

	  overlap.preamble(l) = overlap.preamble(l)+0;

	  ov = (stop_samples - start_samples)/timing.payload_len;
	  overlap.payload(l) = overlap.payload(l)+ov;
	end
        if (stop_samples - start_samples > timing.payload_len)
            fprintf('Alert!!!! stop_samples - start_samples > timing.payload_len\n');
	    exit;
        end

        mui_packet = create_packet(start_samples,stop_samples,compound_channel,...
                                   length(compound_channel),0.5, ...
                                   tx_bit,ep.ths,ep.scrambling,symbol_start,CODE_ID,L, ...
                                   TC,FS_CONT,PACKET_LENGTH,0,PACKET_LENGTH,...
                                   BURST_LENGTH,SYMBOL_LENGTH,...
                                   samples,0);

        fprintf('[%d] Samples of the MUI packet (%d,%d,%d,%d) Payload %d, Packet %d\n',...
                overlap.do(l),...
                start_samples,stop_samples,stop_samples-start_samples,length(mui_packet),...
                timing.payload_len,timing.packet_len);

        rx_data(1:length(mui_packet)) = ...
            rx_data(1:length(mui_packet)) + mui_packet;

      elseif (chunk.mui_delay(l).delay(k) < 0)
        fprintf('MUI packet starts after the data packet\n');
        fprintf('Delay = %12.2f (%.2f)\n',chunk.mui_delay(l).delay(k),...
                1-(-chunk.mui_delay(l).delay(k)/timing.packet_len));

        start_samples = 1;
        stop_samples  = timing.payload_len-abs(chunk.mui_delay(l).delay(k));

	% Stat
	% Do we overlap with the preamble, or the data part. or both?
	if stop_samples < timing.preamble_len
	  fprintf('[%d] Overlap with preamble only\n',overlap.do(l));

	  ov = stop_samples/timing.payload_len;
	  overlap.preamble(l) = overlap.preamble(l)+ov;

	  overlap.payload(l) = overlap.payload(l)+0;
	else
	  fprintf('[%d] Overlap with preamble and data\n',overlap.do(l));

	  overlap.preamble(l) = overlap.preamble(l)+timing.preamble_len/timing.payload_len;

	  ov = (stop_samples - timing.preamble_len)/timing.payload_len;
	  overlap.payload(l) = overlap.payload(l)+ov;
	end

        mui_packet = create_packet(start_samples,stop_samples,compound_channel,...
                                   length(compound_channel),0.5, ...
                                   tx_bit,ep.ths,ep.scrambling,symbol_start,CODE_ID,L, ...
                                   TC,FS_CONT,PACKET_LENGTH,0,PACKET_LENGTH,...
                                   BURST_LENGTH,SYMBOL_LENGTH,...
                                   samples,0);

        fprintf('[%d] Samples of the MUI packet (%d,%d,%d,%d) Payload %d, Packet %d\n',...
                overlap.do(l),...
                start_samples,stop_samples,stop_samples-start_samples,length(mui_packet),...
                timing.payload_len,timing.packet_len);

        % Add the MUI
        rx_data(end-length(mui_packet)+1:end) = ...
            rx_data(end-length(mui_packet)+1:end) + mui_packet;

      end
      % Debug
      fprintf('[%d] Overlap percentage: %.2f %.2f (%.2f)\n',...
              overlap.do(l),...
              overlap.preamble(l),overlap.payload(l),length(mui_packet)/timing.payload_len);

      clear mui_packet;

    end    
  end

  % Add AWGN
  fprintf('add noise, ');
  [noise_0 noise_1] = randn_bl(ALPHA,Ts,Nt,B,s2*(2*B)/FS_CONT);
  %size(rx_data)
  %size([noise_0 noise_1])
  rx_data = rx_data + [noise_0 noise_1];
  clear noise_0 noise_1;

  % Receiver
  %
  fprintf('receive\n');
  demod_bit = ed_decoder(rx_data, ep.ths, ...
                         symbol_start, ...
                         burstlength,samples.burst,samples.symbol, ...
                         channel_mask, channel_mask_unit, ...
                         timing_offset_samples, ...
                         PACKET_LENGTH, FS_CONT, T_INT, N_INT);

  if RSCODE == 1
    % RS decoding
    dec_bit = rs_decoder_154a(demod_bit,ep.interleaver_seed);
  else
    dec_bit = demod_bit;
  end

  % Compute the BER
  num_err(l) = sum(abs(data_bit-dec_bit));
  ber_sim(l) = num_err(l)/DATA_LENGTH;

  % Stat
  overlap.err(l) = num_err(l);
  overlap.ber(l) = ber_sim(l);
end

dt=etime(clock,t0);
fprintf('(dt = %.2f), avg. ber = %.12f\n',dt,mean(ber_sim));
% Save the result in RESULT_FILE
val = [ber_sim' num_err' ...
       overlap.do' overlap.preamble' overlap.payload' ...
       channel_mask_length'];
save(RESULT_FILE,'val','-ascii','-double');


fprintf('Simulation over, bye!\n');