Yujen Ku / Object Detection in Urban Scenarios using Vehicular Communication Signals

Download as
Browse Code ยป

Commit 52917a451fd6a529fabf445636acb7bcb499c0af

Authored by Yujen Ku
1 parent 727ccf17e5
Exists in master

modification of figure code

Showing 1 changed file with 3 additions and 4 deletions Inline Diff

mmwave_objectdetection.m View file @ 52917a4
clear; 1 1 clear;
alpha_est_list = []; 2 2 alpha_est_list = [];
3 3
raw_rx_dec_record = 0; 4 4 raw_rx_dec_record = 0;
%% Params: 5 5 %% Params:
6 6
%control param 7 7 %control param
Detail_plot = 0; 8 8 Detail_plot = 0;
showdetail = 0; 9 9 showdetail = 0;
DETECTION_OFFSET = 50; % to add packet detection error 10 10 DETECTION_OFFSET = 50; % to add packet detection error
INTERPOLATE = 0; 11 11 INTERPOLATE = 0;
12 12
%snr settings (run estimation under different snr values) 13 13 %snr settings (run estimation under different snr values)
snr_value = [0]; %snr_value = [10,20,50,100,200,500,1000]; 14 14 snr_value = [100]; %snr_value = [10,20,50,100,200,500,1000];
15 15
%how many cycles of packets you are going to measure(each cycle: 2 nulling estimation packets andd 5 object detecting packets) 16 16 %how many cycles of packets you are going to measure(each cycle: 2 nulling estimation packets andd 5 object detecting packets)
total_num_est = 100; 17 17 total_num_est = 100;
18 18
%how many packet used in one measurement. (not including 2 packets for nulling) 19 19 %how many packet used in one measurement. (not including 2 packets for nulling)
num_avgpkt = 5; 20 20 num_avgpkt = 5;
21 21
%variables for estimated parameters 22 22 %variables for estimated parameters
obj_est_gridx_record = zeros(length(snr_value),total_num_est); 23 23 obj_est_gridx_record = zeros(length(snr_value),total_num_est);
obj_est_gridy_record = zeros(length(snr_value),total_num_est); 24 24 obj_est_gridy_record = zeros(length(snr_value),total_num_est);
ant_est_veloc_record = zeros(length(snr_value),total_num_est); 25 25 ant_est_veloc_record = zeros(length(snr_value),total_num_est);
26 26
27 27
for snr_i = 1:1:length(snr_value) 28 28 for snr_i = 1:1:length(snr_value)
snr = snr_value(snr_i); 29 29 snr = snr_value(snr_i);
disp([ 'creating case snr = ' num2str(snr)]); 30 30 disp([ 'creating case snr = ' num2str(snr)]);
%environment param 31 31 %environment param
% | || 32 32 % | ||
% (0,0)+-------------------- <---v_obj -----o(x_obj) 33 33 % (0,0)+-------------------- <---v_obj -----o(x_obj)
% | || 34 34 % | ||
% |< x_wall >|| 35 35 % |< x_wall >||
% | || 36 36 % | ||
% | || 37 37 % | ||
% | || 38 38 % | ||
% ^ || 39 39 % ^ ||
% | || 40 40 % | ||
%(y_ant) car(ant) || 41 41 %(y_ant) car(ant) ||
42 42
%setting the initial position, velocity of object and antenna 43 43 %setting the initial position, velocity of object and antenna
x_obj = 3; %initial position of moving object(m) 44 44 x_obj = 3; %initial position of moving object(m)
x_ori_obj = x_obj; 45 45 x_ori_obj = x_obj;
y_ant = -10; %initial position of antenna(m) 46 46 y_ant = -10; %initial position of antenna(m)
v_ant = 10; %velocity of antenna(m/s) 47 47 v_ant = 10; %velocity of antenna(m/s)
v_obj = -1.5; %velocity of moving object(m/s) 48 48 v_obj = -1.5; %velocity of moving object(m/s)
x_wall = 1; %horizontal distance of wall to y axis (m) 49 49 x_wall = 1; %horizontal distance of wall to y axis (m)
thick_wall = 0.01; %2*thickness of wall (m) 50 50 thick_wall = 0.01; %2*thickness of wall (m)
51 51
% __/ Tx 52 52 % __/ Tx
% 53 53 %
% __/ Rx 54 54 % __/ Rx
% 55 55 %
% __/ Tx 56 56 % __/ Tx
57 57
dis_ant = 0.1; %distance between antennas 58 58 dis_ant = 0.1; %distance between antennas
59 59
% Waveform params 60 60 % Waveform params
N_OFDM_SYMS = 500; % Number of OFDM symbols 61 61 N_OFDM_SYMS = 500; % Number of OFDM symbols
MOD_ORDER = 2; % Modulation order in power of 2 (1/2/4/6 = BSPK/QPSK/16-QAM/64-QAM) 62 62 MOD_ORDER = 2; % Modulation order in power of 2 (1/2/4/6 = BSPK/QPSK/16-QAM/64-QAM)
TX_SCALE = 1.0; % Scale for Tx waveform ([0:1]) 63 63 TX_SCALE = 1.0; % Scale for Tx waveform ([0:1])
64 64
% OFDM params 65 65 % OFDM params
SC_IND_PILOTS = [-150,-130,-110,-90,-70,-50,-30,-10,10,30,50,70,90,110,130,150]+256; % Pilot subcarrier indices 66 66 SC_IND_PILOTS = [-150,-130,-110,-90,-70,-50,-30,-10,10,30,50,70,90,110,130,150]+256; % Pilot subcarrier indices
SC_IND_DATA = [-177:-151, -149:-131, -129:-111, -109:-91, -89:-71, -69:-51, -49:-31, -29:-11, -9:-2, 2:9, 11:29, 31:49, 51:69, 71:89, 91:109, 111:129,131:149 , 151:177] + 256; % Data subcarrier indices 67 67 SC_IND_DATA = [-177:-151, -149:-131, -129:-111, -109:-91, -89:-71, -69:-51, -49:-31, -29:-11, -9:-2, 2:9, 11:29, 31:49, 51:69, 71:89, 91:109, 111:129,131:149 , 151:177] + 256; % Data subcarrier indices
N_SC = 512; % Number of subcarriers 68 68 N_SC = 512; % Number of subcarriers
CP_LEN = 128; % Cyclic prefix length 69 69 CP_LEN = 128; % Cyclic prefix length
N_DATA_SYMS = N_OFDM_SYMS * length(SC_IND_DATA); % Number of data symbols (one per data-bearing subcarrier per OFDM symbol) 70 70 N_DATA_SYMS = N_OFDM_SYMS * length(SC_IND_DATA); % Number of data symbols (one per data-bearing subcarrier per OFDM symbol)
71 71
channel_coding = .5; % coding rate 72 72 channel_coding = .5; % coding rate
trellis_end_length = 8; % bits for trellis to end 73 73 trellis_end_length = 8; % bits for trellis to end
74 74
v_c = 299792458; %light speed (m/s) 75 75 v_c = 299792458; %light speed (m/s)
SC_spacing = 5.165625e6; %sub-carrier spacing(Hz) 76 76 SC_spacing = 5.165625e6; %sub-carrier spacing(Hz)
center_freq = 60e9; %center frequency 77 77 center_freq = 60e9; %center frequency
symbol_time = 1/SC_spacing; %duration of 1 OFDM symbol 78 78 symbol_time = 1/SC_spacing; %duration of 1 OFDM symbol
sample_time = symbol_time/N_SC; %duration of 1 sample 79 79 sample_time = symbol_time/N_SC; %duration of 1 sample
80 80
pck_interval = 200e-6; %interval between two packets 81 81 pck_interval = 200e-6; %interval between two packets
sample_interval = ceil(pck_interval/sample_time); %number of samples between two packets 82 82 sample_interval = ceil(pck_interval/sample_time); %number of samples between two packets
83 83
distance_real = zeros(1,total_num_est); %record the real distance of each measurement 84 84 distance_real = zeros(1,total_num_est); %record the real distance of each measurement
distance_esti = zeros(1,total_num_est); %record the estimate distance of each measurement 85 85 distance_esti = zeros(1,total_num_est); %record the estimate distance of each measurement
ant_est_velo = zeros(1,total_num_est); %record the estimate antenna velocity of each measurement 86 86 ant_est_velo = zeros(1,total_num_est); %record the estimate antenna velocity of each measurement
ant_est_velo(1,1) = v_ant; %initial the estimation as the true speed 87 87 ant_est_velo(1,1) = v_ant; %initial the estimation as the true speed
gridx_real = zeros(1,total_num_est); %record the real x position of the object 88 88 gridx_real = zeros(1,total_num_est); %record the real x position of the object
gridy_real = zeros(1,total_num_est); %record the real y position of the object 89 89 gridy_real = zeros(1,total_num_est); %record the real y position of the object
alpha_est_list = zeros(1,total_num_est); %record the estiamte AoA of each measurement 90 90 alpha_est_list = zeros(1,total_num_est); %record the estiamte AoA of each measurement
alpha_real_list = zeros(1,total_num_est); %record the real AoA of each measurement 91 91 alpha_real_list = zeros(1,total_num_est); %record the real AoA of each measurement
est_pktwindow = 30; %do the measurement after every 'est_pktwindow' number of packets 92 92 est_pktwindow = 30; %do the measurement after every 'est_pktwindow' number of packets
num_est = 0; 93 93 num_est = 0;
94 94
obj_mov_ppkt = pck_interval*v_obj; 95 95 obj_mov_ppkt = pck_interval*v_obj;
for t = 1:total_num_est*est_pktwindow 96 96 for t = 1:total_num_est*est_pktwindow
if mod(t,est_pktwindow) == 0 97 97 if mod(t,est_pktwindow) == 0
disp([ ' creating signal t = ' num2str(t)]); 98 98 disp([ ' creating signal t = ' num2str(t)]);
end 99 99 end
100 100
%update the start index of time, sample and position of object 101 101 %update the start index of time, sample and position of object
time_start = (t-1)*pck_interval+1; 102 102 time_start = (t-1)*pck_interval+1;
sample_start = (t-1)*sample_interval+1; 103 103 sample_start = (t-1)*sample_interval+1;
x_obj = x_obj+v_obj*pck_interval; %new x location of object 104 104 x_obj = x_obj+v_obj*pck_interval; %new x location of object
y_ant = y_ant+v_ant*pck_interval; %new y location of antenna 105 105 y_ant = y_ant+v_ant*pck_interval; %new y location of antenna
106 106
%periodically measurement 107 107 %periodically measurement
if mod(t,est_pktwindow)>( 2+num_avgpkt) || mod(t,est_pktwindow) == 0 108 108 if mod(t,est_pktwindow)>( 2+num_avgpkt) || mod(t,est_pktwindow) == 0
%do nothing 109 109 %do nothing
else 110 110 else
111 111
if mod(t,est_pktwindow) == 3 112 112 if mod(t,est_pktwindow) == 3
num_est = num_est+1; 113 113 num_est = num_est+1;
distance_real(1,num_est) = sqrt(x_obj^2+y_ant^2); 114 114 distance_real(1,num_est) = sqrt(x_obj^2+y_ant^2);
end 115 115 end
%% Preamble 116 116 %% Preamble
if mod(t,est_pktwindow) == 1 117 117 if mod(t,est_pktwindow) == 1
H_tx2 = ones(1,N_SC); 118 118 H_tx2 = ones(1,N_SC);
distance_avg = zeros(1,num_avgpkt); 119 119 distance_avg = zeros(1,num_avgpkt);
alpha_avg = zeros(1,num_avgpkt); 120 120 alpha_avg = zeros(1,num_avgpkt);
velo_avg = zeros(1,num_avgpkt); 121 121 velo_avg = zeros(1,num_avgpkt);
num_avgpkt_i = 1; 122 122 num_avgpkt_i = 1;
end 123 123 end
% reference [802.11ad PHY] 124 124 % reference [802.11ad PHY]
125 125
Ga_128_f_1 = [ 1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,1,1,-1,-1,1,1,1,1,-1,1,-1,1,-1,1,1,-1,... 126 126 Ga_128_f_1 = [ 1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,1,1,-1,-1,1,1,1,1,-1,1,-1,1,-1,1,1,-1,...
-1,-1,1,1,1,1,1,1,1,-1,1,-1,-1,1,1,-1,1,1,-1,-1,1,1,1,1,-1,1,-1,1,-1,1,1,-1,... 127 127 -1,-1,1,1,1,1,1,1,1,-1,1,-1,-1,1,1,-1,1,1,-1,-1,1,1,1,1,-1,1,-1,1,-1,1,1,-1,...
1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,1,1,-1,-1,1,1,1,1,-1,1,-1,1,-1,1,1,-1,... 128 128 1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,1,1,-1,-1,1,1,1,1,-1,1,-1,1,-1,1,1,-1,...
1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,-1,-1,1,1,-1,-1,-1,-1,1,-1,1,-1,1,-1,-1,1]; 129 129 1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,-1,-1,1,1,-1,-1,-1,-1,1,-1,1,-1,1,-1,-1,1];
Gb_128_f_1 = [-1,-1, 1, 1, 1, 1, 1, 1, 1,-1, 1,-1,-1, 1, 1,-1,-1,-1,1,1,-1,-1,-1,-1,1,-1,1,-1,1,-1,-1,1,... 130 130 Gb_128_f_1 = [-1,-1, 1, 1, 1, 1, 1, 1, 1,-1, 1,-1,-1, 1, 1,-1,-1,-1,1,1,-1,-1,-1,-1,1,-1,1,-1,1,-1,-1,1,...
1, 1,-1,-1,-1,-1,-1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1,-1,1,1,-1,-1,-1,-1,1,-1,1,-1,1,-1,-1,1,... 131 131 1, 1,-1,-1,-1,-1,-1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1,-1,1,1,-1,-1,-1,-1,1,-1,1,-1,1,-1,-1,1,...
1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,1,1,-1,-1,1,1,1,1,-1,1,-1,1,-1,1,1,-1,... 132 132 1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,1,1,-1,-1,1,1,1,1,-1,1,-1,1,-1,1,1,-1,...
1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,-1,-1,1,1,-1,-1,-1,-1,1,-1,1,-1,1,-1,-1,1]; 133 133 1,1,-1,-1,-1,-1,-1,-1,-1,1,-1,1,1,-1,-1,1,-1,-1,1,1,-1,-1,-1,-1,1,-1,1,-1,1,-1,-1,1];
Ga_128_t = ifft(Ga_128_f_1,128); 134 134 Ga_128_t = ifft(Ga_128_f_1,128);
Ga_128_t_neg = ifft(-1*Ga_128_f_1,128); 135 135 Ga_128_t_neg = ifft(-1*Ga_128_f_1,128);
Gb_128_t = ifft(Gb_128_f_1,128); 136 136 Gb_128_t = ifft(Gb_128_f_1,128);
Gb_128_t_neg = ifft(-1*Gb_128_f_1,128); 137 137 Gb_128_t_neg = ifft(-1*Gb_128_f_1,128);
138 138
Gv_512_f = [-Gb_128_f_1, Ga_128_f_1, -Gb_128_f_1, -Ga_128_f_1]; 139 139 Gv_512_f = [-Gb_128_f_1, Ga_128_f_1, -Gb_128_f_1, -Ga_128_f_1];
Gu_512_f = [-Gb_128_f_1, -Ga_128_f_1, Gb_128_f_1, -Ga_128_f_1]; 140 140 Gu_512_f = [-Gb_128_f_1, -Ga_128_f_1, Gb_128_f_1, -Ga_128_f_1];
141 141
%Nulling for preamble 142 142 %Nulling for preamble
Gv_512_f_2 = H_tx2.*[-Gb_128_f_1, Ga_128_f_1, -Gb_128_f_1, -Ga_128_f_1]; 143 143 Gv_512_f_2 = H_tx2.*[-Gb_128_f_1, Ga_128_f_1, -Gb_128_f_1, -Ga_128_f_1];
Gu_512_f_2 = H_tx2.*[-Gb_128_f_1, -Ga_128_f_1, Gb_128_f_1, -Ga_128_f_1]; 144 144 Gu_512_f_2 = H_tx2.*[-Gb_128_f_1, -Ga_128_f_1, Gb_128_f_1, -Ga_128_f_1];
145 145
%Fourier transform of preamble from freq domain to time domain 146 146 %Fourier transform of preamble from freq domain to time domain
%first antenna 147 147 %first antenna
Gv_512_t = ifft(Gv_512_f,512); 148 148 Gv_512_t = ifft(Gv_512_f,512);
Gu_512_t = ifft(Gu_512_f,512); 149 149 Gu_512_t = ifft(Gu_512_f,512);
%second antenna 150 150 %second antenna
Gv_512_t_2 = ifft(Gv_512_f_2,512); 151 151 Gv_512_t_2 = ifft(Gv_512_f_2,512);
Gu_512_t_2 = ifft(Gu_512_f_2,512); 152 152 Gu_512_t_2 = ifft(Gu_512_f_2,512);
153 153
%produce preamble 154 154 %produce preamble
preamble_1 = [repmat(Ga_128_t, 1, 16) Ga_128_t_neg Gv_512_t Gu_512_t Gb_128_t_neg]; %preamble for Tx1 155 155 preamble_1 = [repmat(Ga_128_t, 1, 16) Ga_128_t_neg Gv_512_t Gu_512_t Gb_128_t_neg]; %preamble for Tx1
preamble_2 = [repmat(Ga_128_t, 1, 16) Ga_128_t_neg Gv_512_t_2 Gu_512_t_2 Gb_128_t_neg]; %preamble for Tx2 156 156 preamble_2 = [repmat(Ga_128_t, 1, 16) Ga_128_t_neg Gv_512_t_2 Gu_512_t_2 Gb_128_t_neg]; %preamble for Tx2
157 157
%% Generate a payload of random integers 158 158 %% Generate a payload of random integers
number_of_bits = (N_DATA_SYMS * MOD_ORDER - 2*trellis_end_length) * channel_coding; 159 159 number_of_bits = (N_DATA_SYMS * MOD_ORDER - 2*trellis_end_length) * channel_coding;
tx_data = randi(2, 1, number_of_bits) - 1; 160 160 tx_data = randi(2, 1, number_of_bits) - 1;
161 161
% Forward Error Correction 162 162 % Forward Error Correction
tx_data = double([tx_data zeros(1,trellis_end_length) ]); % 8 bits padding 163 163 tx_data = double([tx_data zeros(1,trellis_end_length) ]); % 8 bits padding
trel = poly2trellis(7, [171 133]); % Define trellis 164 164 trel = poly2trellis(7, [171 133]); % Define trellis
tx_code = convenc(tx_data,trel); % convultional encoder 165 165 tx_code = convenc(tx_data,trel); % convultional encoder
166 166
% bits to signal space mapping these are your x_k from the class 167 167 % bits to signal space mapping these are your x_k from the class
tx_syms = mapping(tx_code', MOD_ORDER, 1); 168 168 tx_syms = mapping(tx_code', MOD_ORDER, 1);
169 169
if Detail_plot 170 170 if Detail_plot
figure(1); 171 171 figure(1);
scatter(real(tx_syms), imag(tx_syms),'filled'); 172 172 scatter(real(tx_syms), imag(tx_syms),'filled');
title(' Signal Space of transmitted bits'); 173 173 title(' Signal Space of transmitted bits');
xlabel('I'); ylabel('Q'); 174 174 xlabel('I'); ylabel('Q');
title('Tx data for 64QAM data set') 175 175 title('Tx data for 64QAM data set')
end 176 176 end
177 177
% Reshape the symbol vector to a matrix with one column per OFDM symbol, 178 178 % Reshape the symbol vector to a matrix with one column per OFDM symbol,
tx_syms_mat = reshape(tx_syms, length(SC_IND_DATA), N_OFDM_SYMS); 179 179 tx_syms_mat = reshape(tx_syms, length(SC_IND_DATA), N_OFDM_SYMS);
180 180
% Define the pilot tone values as BPSK symbols 181 181 % Define the pilot tone values as BPSK symbols
pilots = [-1, 1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, 1, 1].'; 182 182 pilots = [-1, 1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, 1, 1].';
183 183
% Repeat the pilots across all OFDM symbols 184 184 % Repeat the pilots across all OFDM symbols
pilots_mat = repmat(pilots, 1, N_OFDM_SYMS); 185 185 pilots_mat = repmat(pilots, 1, N_OFDM_SYMS);
186 186
187 187
%% IFFT 188 188 %% IFFT
189 189
% Construct the IFFT input matrix 190 190 % Construct the IFFT input matrix
ifft_in_mat = zeros(N_SC, N_OFDM_SYMS); 191 191 ifft_in_mat = zeros(N_SC, N_OFDM_SYMS);
192 192
% Insert the data and pilot values; other subcarriers will remain at 0 193 193 % Insert the data and pilot values; other subcarriers will remain at 0
ifft_in_mat(SC_IND_DATA, :) = tx_syms_mat; 194 194 ifft_in_mat(SC_IND_DATA, :) = tx_syms_mat;
ifft_in_mat(SC_IND_PILOTS, :) = pilots_mat; 195 195 ifft_in_mat(SC_IND_PILOTS, :) = pilots_mat;
196 196
%Perform the IFFT --> frequency to time translation 197 197 %Perform the IFFT --> frequency to time translation
tx_payload_mat = ifft(ifft_in_mat, N_SC, 1); %payload for Tx1 198 198 tx_payload_mat = ifft(ifft_in_mat, N_SC, 1); %payload for Tx1
tx_payload_mat_2 = ifft(ifft_in_mat.*repmat(H_tx2',1,N_OFDM_SYMS), N_SC, 1); %payload for Tx2 199 199 tx_payload_mat_2 = ifft(ifft_in_mat.*repmat(H_tx2',1,N_OFDM_SYMS), N_SC, 1); %payload for Tx2
% Insert the cyclic prefix 200 200 % Insert the cyclic prefix
if(CP_LEN > 0) 201 201 if(CP_LEN > 0)
tx_cp = tx_payload_mat((end-CP_LEN+1 : end), :); 202 202 tx_cp = tx_payload_mat((end-CP_LEN+1 : end), :);
tx_payload_mat = [tx_cp; tx_payload_mat]; 203 203 tx_payload_mat = [tx_cp; tx_payload_mat];
tx_cp_2 = tx_payload_mat_2((end-CP_LEN+1 : end), :); 204 204 tx_cp_2 = tx_payload_mat_2((end-CP_LEN+1 : end), :);
tx_payload_mat_2 = [tx_cp_2; tx_payload_mat_2]; 205 205 tx_payload_mat_2 = [tx_cp_2; tx_payload_mat_2];
end 206 206 end
207 207
% Reshape to a vector 208 208 % Reshape to a vector
tx_payload_vec = reshape(tx_payload_mat, 1, numel(tx_payload_mat)); 209 209 tx_payload_vec = reshape(tx_payload_mat, 1, numel(tx_payload_mat));
tx_payload_vec_2 = reshape(tx_payload_mat_2, 1, numel(tx_payload_mat_2)); 210 210 tx_payload_vec_2 = reshape(tx_payload_mat_2, 1, numel(tx_payload_mat_2));
211 211
% Construct the full time-domain OFDM waveform 212 212 % Construct the full time-domain OFDM waveform
tx_vec = [preamble_1 tx_payload_vec]; 213 213 tx_vec = [preamble_1 tx_payload_vec];
tx_vec_2 = [preamble_2 tx_payload_vec_2]; 214 214 tx_vec_2 = [preamble_2 tx_payload_vec_2];
215 215
%% Interpolate, 216 216 %% Interpolate,
if INTERPOLATE 217 217 if INTERPOLATE
% Interpolation filter basically implements the DAC before transmission 218 218 % Interpolation filter basically implements the DAC before transmission
% On the receiver's end decimation is performed to implement the ADC 219 219 % On the receiver's end decimation is performed to implement the ADC
220 220
% Define a half-band 2x interpolation filter response 221 221 % Define a half-band 2x interpolation filter response
interp_filt2 = zeros(1,43); 222 222 interp_filt2 = zeros(1,43);
interp_filt2([1 3 5 7 9 11 13 15 17 19 21]) = [12 -32 72 -140 252 -422 682 -1086 1778 -3284 10364]; 223 223 interp_filt2([1 3 5 7 9 11 13 15 17 19 21]) = [12 -32 72 -140 252 -422 682 -1086 1778 -3284 10364];
interp_filt2([23 25 27 29 31 33 35 37 39 41 43]) = interp_filt2(fliplr([1 3 5 7 9 11 13 15 17 19 21])); 224 224 interp_filt2([23 25 27 29 31 33 35 37 39 41 43]) = interp_filt2(fliplr([1 3 5 7 9 11 13 15 17 19 21]));
interp_filt2(22) = 16384; 225 225 interp_filt2(22) = 16384;
interp_filt2 = interp_filt2./max(abs(interp_filt2)); 226 226 interp_filt2 = interp_filt2./max(abs(interp_filt2));
227 227
% Pad with zeros for transmission to deal with delay through the interpolation filter 228 228 % Pad with zeros for transmission to deal with delay through the interpolation filter
tx_vec_padded = [tx_vec, zeros(1, ceil(length(interp_filt2)/2))]; 229 229 tx_vec_padded = [tx_vec, zeros(1, ceil(length(interp_filt2)/2))];
tx_vec_2x = zeros(1, 2*numel(tx_vec_padded)); 230 230 tx_vec_2x = zeros(1, 2*numel(tx_vec_padded));
tx_vec_2x(1:2:end) = tx_vec_padded; 231 231 tx_vec_2x(1:2:end) = tx_vec_padded;
tx_vec_air = filter(interp_filt2, 1, tx_vec_2x); 232 232 tx_vec_air = filter(interp_filt2, 1, tx_vec_2x);
233 233
if Detail_plot 234 234 if Detail_plot
figure(2); 235 235 figure(2);
plot(abs(tx_vec_2x)); 236 236 plot(abs(tx_vec_2x));
hold on; 237 237 hold on;
plot(abs(tx_vec_air(22:end))); 238 238 plot(abs(tx_vec_air(22:end)));
xlim([20,50]); 239 239 xlim([20,50]);
title('Interpolation visualized'); 240 240 title('Interpolation visualized');
xlabel('time'); ylabel('amplitude'); 241 241 xlabel('time'); ylabel('amplitude');
legend('y = pre filtering','y = post filtering') 242 242 legend('y = pre filtering','y = post filtering')
end 243 243 end
% Scale the Tx vector to +/- 1, becasue ADC and DAC take samples input from 244 244 % Scale the Tx vector to +/- 1, becasue ADC and DAC take samples input from
% 1 to -1 245 245 % 1 to -1
tx_vec_air = TX_SCALE .* tx_vec_air ./ max(abs(tx_vec_air)); 246 246 tx_vec_air = TX_SCALE .* tx_vec_air ./ max(abs(tx_vec_air));
247 247
if Detail_plot 248 248 if Detail_plot
figure(3); 249 249 figure(3);
plot(db(abs(fftshift(fft(tx_vec_air))))); 250 250 plot(db(abs(fftshift(fft(tx_vec_air)))));
%plot(db(abs(fftshift(fft(tx_vec_2x))))); 251 251 %plot(db(abs(fftshift(fft(tx_vec_2x)))));
%xlim([20000,60000]); ylim([0,65]); 252 252 %xlim([20000,60000]); ylim([0,65]);
% in this plot, why do see four peaks? 253 253 % in this plot, why do see four peaks?
end 254 254 end
else 255 255 else
tx_vec_air = tx_vec; 256 256 tx_vec_air = tx_vec;
tx_vec_air_2 = tx_vec_2; 257 257 tx_vec_air_2 = tx_vec_2;
end 258 258 end
%% This part of code is for simulating the wireless channel. 259 259 %% This part of code is for simulating the wireless channel.
260 260
%send to 2 transmitter: 261 261 %send to 2 transmitter:
tx1_vec_air = tx_vec_air; 262 262 tx1_vec_air = tx_vec_air;
tx2_vec_air = tx_vec_air_2; 263 263 tx2_vec_air = tx_vec_air_2;
264 264
% AWGN: 265 265 % AWGN:
TRIGGER_OFFSET_TOL_NS = 3000; % Trigger time offset toleration between Tx and Rx that can be accomodated 266 266 TRIGGER_OFFSET_TOL_NS = 3000; % Trigger time offset toleration between Tx and Rx that can be accomodated
SAMP_FREQ = 1/sample_time; % Sampling frequency 267 267 SAMP_FREQ = 1/sample_time; % Sampling frequency
268 268
rx_vec_air_tx1 = [tx1_vec_air, zeros(1,ceil((TRIGGER_OFFSET_TOL_NS*1e-9) / (1/SAMP_FREQ)))]; 269 269 rx_vec_air_tx1 = [tx1_vec_air, zeros(1,ceil((TRIGGER_OFFSET_TOL_NS*1e-9) / (1/SAMP_FREQ)))];
rx_vec_air_tx2 = [tx2_vec_air, zeros(1,ceil((TRIGGER_OFFSET_TOL_NS*1e-9) / (1/SAMP_FREQ)))]; 270 270 rx_vec_air_tx2 = [tx2_vec_air, zeros(1,ceil((TRIGGER_OFFSET_TOL_NS*1e-9) / (1/SAMP_FREQ)))];
271 271
%add noise 272 272 %add noise
noise_power = var(rx_vec_air_tx1) * 10 ^(-snr/20); 273 273 noise_power = var(rx_vec_air_tx1) * 10 ^(-snr/20);
rx_vec_air_tx1 = rx_vec_air_tx1 + noise_power*complex(randn(1,length(rx_vec_air_tx1)), randn(1,length(rx_vec_air_tx1))); 274 274 rx_vec_air_tx1 = rx_vec_air_tx1 + noise_power*complex(randn(1,length(rx_vec_air_tx1)), randn(1,length(rx_vec_air_tx1)));
275 275
noise_power = var(rx_vec_air_tx2) * 10 ^(-snr/20); 276 276 noise_power = var(rx_vec_air_tx2) * 10 ^(-snr/20);
rx_vec_air_tx2 = rx_vec_air_tx2 + noise_power*complex(randn(1,length(rx_vec_air_tx2)), randn(1,length(rx_vec_air_tx2))); 277 277 rx_vec_air_tx2 = rx_vec_air_tx2 + noise_power*complex(randn(1,length(rx_vec_air_tx2)), randn(1,length(rx_vec_air_tx2)));
278 278
% Decimate %DAC in receiver 279 279 % Decimate %DAC in receiver
if INTERPOLATE 280 280 if INTERPOLATE
raw_rx_dec_tx1 = filter(interp_filt2, 1, rx_vec_air_tx1); 281 281 raw_rx_dec_tx1 = filter(interp_filt2, 1, rx_vec_air_tx1);
raw_rx_dec_tx1 = [zeros(1,DETECTION_OFFSET) raw_rx_dec_tx1(1:2:end)]; 282 282 raw_rx_dec_tx1 = [zeros(1,DETECTION_OFFSET) raw_rx_dec_tx1(1:2:end)];
283 283
raw_rx_dec_tx2 = filter(interp_filt2, 1, rx_vec_air_tx2); 284 284 raw_rx_dec_tx2 = filter(interp_filt2, 1, rx_vec_air_tx2);
raw_rx_dec_tx2 = [zeros(1,DETECTION_OFFSET) raw_rx_dec_tx2(1:2:end)]; 285 285 raw_rx_dec_tx2 = [zeros(1,DETECTION_OFFSET) raw_rx_dec_tx2(1:2:end)];
else 286 286 else
raw_rx_dec_tx1 = rx_vec_air_tx1; 287 287 raw_rx_dec_tx1 = rx_vec_air_tx1;
raw_rx_dec_tx2 = rx_vec_air_tx2; 288 288 raw_rx_dec_tx2 = rx_vec_air_tx2;
end 289 289 end
290 290
%create four paths and their delay, phase change, doppler 291 291 %create four paths and their delay, phase change, doppler
%calculating distance 292 292 %calculating distance
tempa_1 = (-2*dis_ant-sqrt( 4*dis_ant^2 - 4*(1-v_c^2/v_ant^2)*(4*x_wall^2+dis_ant^2) ))/2/(1-v_c^2/v_ant^2); %precise 293 293 tempa_1 = (-2*dis_ant-sqrt( 4*dis_ant^2 - 4*(1-v_c^2/v_ant^2)*(4*x_wall^2+dis_ant^2) ))/2/(1-v_c^2/v_ant^2); %precise
tempa_2 = (-2*dis_ant+sqrt( 4*dis_ant^2 - 4*(1-v_c^2/v_ant^2)*(4*x_wall^2+dis_ant^2) ))/2/(1-v_c^2/v_ant^2); %precise 294 294 tempa_2 = (-2*dis_ant+sqrt( 4*dis_ant^2 - 4*(1-v_c^2/v_ant^2)*(4*x_wall^2+dis_ant^2) ))/2/(1-v_c^2/v_ant^2); %precise
distance_tx1wal = tempa_1*v_c/v_ant; 295 295 distance_tx1wal = tempa_1*v_c/v_ant;
distance_tx1obj = sqrt(x_obj^2+(y_ant-dis_ant)^2)+sqrt(x_obj^2+(y_ant-tempa_2)^2); 296 296 distance_tx1obj = sqrt(x_obj^2+(y_ant-dis_ant)^2)+sqrt(x_obj^2+(y_ant-tempa_2)^2);
distance_tx2wal = -tempa_2*v_c/v_ant; 297 297 distance_tx2wal = -tempa_2*v_c/v_ant;
distance_tx2obj = sqrt(x_obj^2+(y_ant+dis_ant)^2)+sqrt(x_obj^2+(y_ant-tempa_2)^2); 298 298 distance_tx2obj = sqrt(x_obj^2+(y_ant+dis_ant)^2)+sqrt(x_obj^2+(y_ant-tempa_2)^2);
299 299
%calculating correspoding delay 300 300 %calculating correspoding delay
Delay_tx1wal = distance_tx1wal/v_c; 301 301 Delay_tx1wal = distance_tx1wal/v_c;
Delay_tx1obj = distance_tx1obj/v_c; 302 302 Delay_tx1obj = distance_tx1obj/v_c;
Delay_tx2wal = distance_tx2wal/v_c; 303 303 Delay_tx2wal = distance_tx2wal/v_c;
Delay_tx2obj = distance_tx2obj/v_c; 304 304 Delay_tx2obj = distance_tx2obj/v_c;
305 305
%calculating channel condition(phase shift) 306 306 %calculating channel condition(phase shift)
Phase_tx1wal = exp(1i*2*pi*distance_tx1wal*center_freq/v_c); 307 307 Phase_tx1wal = exp(1i*2*pi*distance_tx1wal*center_freq/v_c);
Phase_tx1obj = exp(1i*2*pi*distance_tx1obj*center_freq/v_c); 308 308 Phase_tx1obj = exp(1i*2*pi*distance_tx1obj*center_freq/v_c);
Phase_tx2wal = exp(1i*2*pi*distance_tx2wal*center_freq/v_c); 309 309 Phase_tx2wal = exp(1i*2*pi*distance_tx2wal*center_freq/v_c);
Phase_tx2obj = exp(1i*2*pi*distance_tx2obj*center_freq/v_c); 310 310 Phase_tx2obj = exp(1i*2*pi*distance_tx2obj*center_freq/v_c);
311 311
%attenuation according to distance 312 312 %attenuation according to distance
Atten_tx1wal = 68.0630 + 20*log10(distance_tx1wal);%dB 313 313 Atten_tx1wal = 68.0630 + 20*log10(distance_tx1wal);%dB
Atten_tx1obj_wal = 400*( thick_wall*abs( y_ant/x_obj ) ) ;%dB 314 314 Atten_tx1obj_wal = 400*( thick_wall*abs( y_ant/x_obj ) ) ;%dB
Atten_tx1obj_air = 68.0630 + 20*log10(distance_tx1obj-thick_wall*abs(y_ant/x_obj));%dB 315 315 Atten_tx1obj_air = 68.0630 + 20*log10(distance_tx1obj-thick_wall*abs(y_ant/x_obj));%dB
Atten_tx2wal = 68.0630 + 20*log10(distance_tx2wal);%dB 316 316 Atten_tx2wal = 68.0630 + 20*log10(distance_tx2wal);%dB
Atten_tx2obj_wal = 400*( thick_wall*abs(y_ant/x_obj) ) ;%dB 317 317 Atten_tx2obj_wal = 400*( thick_wall*abs(y_ant/x_obj) ) ;%dB
Atten_tx2obj_air = 68.0630 + 20*log10(distance_tx2obj-thick_wall*abs(y_ant/x_obj));%dB 318 318 Atten_tx2obj_air = 68.0630 + 20*log10(distance_tx2obj-thick_wall*abs(y_ant/x_obj));%dB
319 319
%Doppler effect(Frequece shift) 320 320 %Doppler effect(Frequece shift)
theda_velo = atan(abs(v_ant/v_obj)); 321 321 theda_velo = atan(abs(v_ant/v_obj));
theda_grid = atan(abs(y_ant/x_obj)); 322 322 theda_grid = atan(abs(y_ant/x_obj));
theda_dopr = abs(theda_velo-theda_grid); 323 323 theda_dopr = abs(theda_velo-theda_grid);
doppler_FO = center_freq* sqrt( v_ant^2 + v_obj^2 )* cos(theda_dopr)/v_c ; 324 324 doppler_FO = center_freq* sqrt( v_ant^2 + v_obj^2 )* cos(theda_dopr)/v_c ;
325 325
%create 4 signals for each path 326 326 %create 4 signals for each path
raw_rx_dec_tx1wal = raw_rx_dec_tx1*Phase_tx1wal/sqrt( 10^(Atten_tx1wal/10) ); 327 327 raw_rx_dec_tx1wal = raw_rx_dec_tx1*Phase_tx1wal/sqrt( 10^(Atten_tx1wal/10) );
raw_rx_dec_tx1obj = raw_rx_dec_tx1*Phase_tx1obj/sqrt( 10^((Atten_tx1obj_wal+Atten_tx1obj_air)/10) ).* exp(-1i*2*pi*doppler_FO*sample_time*[0:length(raw_rx_dec_tx1)-1]); 328 328 raw_rx_dec_tx1obj = raw_rx_dec_tx1*Phase_tx1obj/sqrt( 10^((Atten_tx1obj_wal+Atten_tx1obj_air)/10) ).* exp(-1i*2*pi*doppler_FO*sample_time*[0:length(raw_rx_dec_tx1)-1]);
raw_rx_dec_tx2wal = raw_rx_dec_tx2*Phase_tx2wal/sqrt( 10^(Atten_tx2wal/10) ); 329 329 raw_rx_dec_tx2wal = raw_rx_dec_tx2*Phase_tx2wal/sqrt( 10^(Atten_tx2wal/10) );
raw_rx_dec_tx2obj = raw_rx_dec_tx2*Phase_tx2obj/sqrt( 10^((Atten_tx2obj_wal+Atten_tx2obj_air)/10) ).* exp(-1i*2*pi*doppler_FO*sample_time*[0:length(raw_rx_dec_tx2)-1]); 330 330 raw_rx_dec_tx2obj = raw_rx_dec_tx2*Phase_tx2obj/sqrt( 10^((Atten_tx2obj_wal+Atten_tx2obj_air)/10) ).* exp(-1i*2*pi*doppler_FO*sample_time*[0:length(raw_rx_dec_tx2)-1]);
331 331
%conbine all 4 reflected signal in reciver 332 332 %conbine all 4 reflected signal in reciver
raw_rx_dec = 0; 333 333 raw_rx_dec = 0;
if mod(t,est_pktwindow) ~= 2 %at step 2, there is no signal from TX1 334 334 if mod(t,est_pktwindow) ~= 2 %at step 2, there is no signal from TX1
if showdetail 335 335 if showdetail
fprintf('t = %d , Wall Reflc-signal from Tx1 with ', t ); 336 336 fprintf('t = %d , Wall Reflc-signal from Tx1 with ', t );
end 337 337 end
raw_rx_dec = form_raw_rx_dec(1,raw_rx_dec,Delay_tx1wal,sample_time,raw_rx_dec_tx1wal,center_freq,showdetail); % refleted signal from tx 1 to wall to Rx 338 338 raw_rx_dec = form_raw_rx_dec(1,raw_rx_dec,Delay_tx1wal,sample_time,raw_rx_dec_tx1wal,center_freq,showdetail); % refleted signal from tx 1 to wall to Rx
if showdetail 339 339 if showdetail
fprintf('t = %d , Objt Reflc-signal from Tx1 with ', t ); 340 340 fprintf('t = %d , Objt Reflc-signal from Tx1 with ', t );
end 341 341 end
raw_rx_dec = form_raw_rx_dec(1,raw_rx_dec,Delay_tx1obj,sample_time,raw_rx_dec_tx1obj,center_freq,showdetail); % refleted signal from tx 1 to object to Rx 342 342 raw_rx_dec = form_raw_rx_dec(1,raw_rx_dec,Delay_tx1obj,sample_time,raw_rx_dec_tx1obj,center_freq,showdetail); % refleted signal from tx 1 to object to Rx
end 343 343 end
if mod(t,est_pktwindow) ~= 1 %at step 1, there is no signal from TX2 344 344 if mod(t,est_pktwindow) ~= 1 %at step 1, there is no signal from TX2
if showdetail 345 345 if showdetail
fprintf('t = %d , Wall Reflc-signal from Tx2 with ', t ); 346 346 fprintf('t = %d , Wall Reflc-signal from Tx2 with ', t );
end 347 347 end
raw_rx_dec = form_raw_rx_dec(1,raw_rx_dec,Delay_tx2wal,sample_time,raw_rx_dec_tx2wal,center_freq,showdetail); % refleted signal from tx 2 to wall to Rx 348 348 raw_rx_dec = form_raw_rx_dec(1,raw_rx_dec,Delay_tx2wal,sample_time,raw_rx_dec_tx2wal,center_freq,showdetail); % refleted signal from tx 2 to wall to Rx
if showdetail 349 349 if showdetail
fprintf('t = %d , Objt Reflc-signal from Tx2 with ', t ); 350 350 fprintf('t = %d , Objt Reflc-signal from Tx2 with ', t );
end 351 351 end
raw_rx_dec = form_raw_rx_dec(1,raw_rx_dec,Delay_tx2obj,sample_time,raw_rx_dec_tx2obj,center_freq,showdetail); 352 352 raw_rx_dec = form_raw_rx_dec(1,raw_rx_dec,Delay_tx2obj,sample_time,raw_rx_dec_tx2obj,center_freq,showdetail);
end 353 353 end
354 354
%% Reveiver Ends 355 355 %% Reveiver Ends
if mod(t,est_pktwindow) == 1 || mod(t,est_pktwindow) == 2 %step 1 and step 2 356 356 if mod(t,est_pktwindow) == 1 || mod(t,est_pktwindow) == 2 %step 1 and step 2
%packet detection 357 357 %packet detection
lts_corr = abs(conv(conj(fliplr(Ga_128_t)), sign(raw_rx_dec))); 358 358 lts_corr = abs(conv(conj(fliplr(Ga_128_t)), sign(raw_rx_dec)));
lts_corr = lts_corr(32:end-32); 359 359 lts_corr = lts_corr(32:end-32);
LTS_CORR_THRESH=.8; 360 360 LTS_CORR_THRESH=.8;
lts_peaks = find(lts_corr > LTS_CORR_THRESH*max(lts_corr)); 361 361 lts_peaks = find(lts_corr > LTS_CORR_THRESH*max(lts_corr));
[LTS1, LTS2] = meshgrid(lts_peaks,lts_peaks); 362 362 [LTS1, LTS2] = meshgrid(lts_peaks,lts_peaks);
[lts_last_peak_index,y] = find(LTS2-LTS1 == length(Ga_128_t)); 363 363 [lts_last_peak_index,y] = find(LTS2-LTS1 == length(Ga_128_t));
364 364
channel_training_ind = lts_peaks(max(lts_last_peak_index)) + 128/2; 365 365 channel_training_ind = lts_peaks(max(lts_last_peak_index)) + 128/2;
payload_ind = lts_peaks(max(lts_last_peak_index)) + 128/2 + 9*128; 366 366 payload_ind = lts_peaks(max(lts_last_peak_index)) + 128/2 + 9*128;
367 367
%CFO 368 368 %CFO
Ga_128_t_start_ind = min(max( lts_peaks(max(lts_last_peak_index)) - 16.5*128,1),length(raw_rx_dec)-16.5*128); 369 369 Ga_128_t_start_ind = min(max( lts_peaks(max(lts_last_peak_index)) - 16.5*128,1),length(raw_rx_dec)-16.5*128);
rx_Ga_128_t = raw_rx_dec(Ga_128_t_start_ind:Ga_128_t_start_ind+16.5*128-1); 370 370 rx_Ga_128_t = raw_rx_dec(Ga_128_t_start_ind:Ga_128_t_start_ind+16.5*128-1);
rx_cfo_est = zeros(1,15); 371 371 rx_cfo_est = zeros(1,15);
for iGa = 1:15 372 372 for iGa = 1:15
rx_Ga_128_1 = rx_Ga_128_t(1+(iGa-1)*128:iGa*128); 373 373 rx_Ga_128_1 = rx_Ga_128_t(1+(iGa-1)*128:iGa*128);
rx_Ga_128_2 = rx_Ga_128_t(1+iGa*128:(iGa+1)*128); 374 374 rx_Ga_128_2 = rx_Ga_128_t(1+iGa*128:(iGa+1)*128);
rx_cfo_est(1,iGa) = mean(unwrap(angle(rx_Ga_128_2 .* conj(rx_Ga_128_1)))); 375 375 rx_cfo_est(1,iGa) = mean(unwrap(angle(rx_Ga_128_2 .* conj(rx_Ga_128_1))));
rx_cfo_est(1,iGa) = rx_cfo_est(1,iGa)/(2*pi*128); 376 376 rx_cfo_est(1,iGa) = rx_cfo_est(1,iGa)/(2*pi*128);
end 377 377 end
rx_cfo_est_final = mean(rx_cfo_est); 378 378 rx_cfo_est_final = mean(rx_cfo_est);
379 379
rx_cfo_corr_t = exp(-1i*2*pi*rx_cfo_est_final*[0:length(raw_rx_dec)-1]); 380 380 rx_cfo_corr_t = exp(-1i*2*pi*rx_cfo_est_final*[0:length(raw_rx_dec)-1]);
rx_dec_cfo_corr = raw_rx_dec .* rx_cfo_corr_t; 381 381 rx_dec_cfo_corr = raw_rx_dec .* rx_cfo_corr_t;
382 382
% channel estimation: 383 383 % channel estimation:
Gv512_ind_start = channel_training_ind ; 384 384 Gv512_ind_start = channel_training_ind ;
Gv512_ind_end = channel_training_ind + 512 - 1; 385 385 Gv512_ind_end = channel_training_ind + 512 - 1;
Gu512_ind_start = channel_training_ind + 512; 386 386 Gu512_ind_start = channel_training_ind + 512;
Gu512_ind_end = channel_training_ind + 512*2 - 1; 387 387 Gu512_ind_end = channel_training_ind + 512*2 - 1;
388 388
rx_Gv512 = raw_rx_dec(Gv512_ind_start:Gv512_ind_end); 389 389 rx_Gv512 = raw_rx_dec(Gv512_ind_start:Gv512_ind_end);
rx_Gu512 = raw_rx_dec(Gu512_ind_start:Gu512_ind_end); 390 390 rx_Gu512 = raw_rx_dec(Gu512_ind_start:Gu512_ind_end);
391 391
rx_Gv512_f = fft(rx_Gv512); 392 392 rx_Gv512_f = fft(rx_Gv512);
rx_Gu512_f = fft(rx_Gu512); 393 393 rx_Gu512_f = fft(rx_Gu512);
394 394
H_vest = rx_Gv512_f./Gv_512_f; 395 395 H_vest = rx_Gv512_f./Gv_512_f;
H_uest = rx_Gu512_f./Gu_512_f; 396 396 H_uest = rx_Gu512_f./Gu_512_f;
H_est = ( H_vest+H_uest )/2; 397 397 H_est = ( H_vest+H_uest )/2;
if mod(t,est_pktwindow) == 1 398 398 if mod(t,est_pktwindow) == 1
H_1 = H_est; 399 399 H_1 = H_est;
else 400 400 else
H_tx2 = - H_1./H_est; 401 401 H_tx2 = - H_1./H_est;
end 402 402 end
lts_corr = abs(conv(conj(fliplr(Gu_512_t)), sign(raw_rx_dec))); 403 403 lts_corr = abs(conv(conj(fliplr(Gu_512_t)), sign(raw_rx_dec)));
lts_corr = lts_corr(32:end-32); 404 404 lts_corr = lts_corr(32:end-32);
LTS_CORR_THRESH=.8; 405 405 LTS_CORR_THRESH=.8;
lts_peaks = find(lts_corr > LTS_CORR_THRESH*max(lts_corr)); 406 406 lts_peaks = find(lts_corr > LTS_CORR_THRESH*max(lts_corr));
407 407
else %mod(t,num_avgpkt) == 0 || 3~6: step 3 408 408 else %mod(t,num_avgpkt) == 0 || 3~6: step 3
%% start extract object information 409 409 %% start extract object information
%packet detection 410 410 %packet detection
corr_target = [ Gv_512_t, Gu_512_t]; 411 411 corr_target = [ Gv_512_t, Gu_512_t];
lts_corr = abs(conv(conj(fliplr(corr_target)), sign(raw_rx_dec))); 412 412 lts_corr = abs(conv(conj(fliplr(corr_target)), sign(raw_rx_dec)));
lts_corr = lts_corr(32:end-32); 413 413 lts_corr = lts_corr(32:end-32);
LTS_CORR_THRESH=.8; 414 414 LTS_CORR_THRESH=.8;
[~,lts_peaks] = max(lts_corr); 415 415 [~,lts_peaks] = max(lts_corr);
416 416
%ToF estimation 417 417 %ToF estimation
D_3 = lts_peaks-128*17-1024+32 ; 418 418 D_3 = lts_peaks-128*17-1024+32 ;
distance_avg(num_avgpkt_i) = D_3*sample_time*v_c/2; 419 419 distance_avg(num_avgpkt_i) = D_3*sample_time*v_c/2;
%average smoothing 420 420 %average smoothing
if num_avgpkt_i == 5 421 421 if num_avgpkt_i == 5
distance_esti(1,num_est) = mean(distance_avg) ; 422 422 distance_esti(1,num_est) = mean(distance_avg) ;
end 423 423 end
if showdetail 424 424 if showdetail
disp(['OBJT delay estimation = ' num2str(D_3)]); 425 425 disp(['OBJT delay estimation = ' num2str(D_3)]);
end 426 426 end
427 427
%CFO estimation: 428 428 %CFO estimation:
Ga_128_t_start_ind = min ( max(lts_peaks-1024 -128*17,1), length(raw_rx_dec) - 16*128 +1) ; 429 429 Ga_128_t_start_ind = min ( max(lts_peaks-1024 -128*17,1), length(raw_rx_dec) - 16*128 +1) ;
rx_Ga_128_t = raw_rx_dec(Ga_128_t_start_ind:Ga_128_t_start_ind+16*128-1); 430 430 rx_Ga_128_t = raw_rx_dec(Ga_128_t_start_ind:Ga_128_t_start_ind+16*128-1);
rx_cfo_est = zeros(1,15); 431 431 rx_cfo_est = zeros(1,15);
for iGa = 1:15 432 432 for iGa = 1:15
rx_Ga_128_1 = rx_Ga_128_t(1+(iGa-1)*128:iGa*128); 433 433 rx_Ga_128_1 = rx_Ga_128_t(1+(iGa-1)*128:iGa*128);
rx_Ga_128_2 = rx_Ga_128_t(1+iGa*128:(iGa+1)*128); 434 434 rx_Ga_128_2 = rx_Ga_128_t(1+iGa*128:(iGa+1)*128);
rx_cfo_est(1,iGa) = mean(unwrap(angle(rx_Ga_128_2 .* conj(rx_Ga_128_1)))); 435 435 rx_cfo_est(1,iGa) = mean(unwrap(angle(rx_Ga_128_2 .* conj(rx_Ga_128_1))));
rx_cfo_est(1,iGa) = rx_cfo_est(1,iGa)/(2*pi*128); 436 436 rx_cfo_est(1,iGa) = rx_cfo_est(1,iGa)/(2*pi*128);
end 437 437 end
rx_cfo_est_final = mean(rx_cfo_est); 438 438 rx_cfo_est_final = mean(rx_cfo_est);
%antenna speed estimation 439 439 %antenna speed estimation
CFO_v_estimate = -rx_cfo_est_final/sample_time*v_c/center_freq; 440 440 CFO_v_estimate = -rx_cfo_est_final/sample_time*v_c/center_freq;
if CFO_v_estimate-v_ant>100 || CFO_v_estimate<0 441 441 if CFO_v_estimate-v_ant>100 || CFO_v_estimate<0
velo_avg(num_avgpkt_i) = ant_est_velo(1,max( num_est-1, 1)); 442 442 velo_avg(num_avgpkt_i) = ant_est_velo(1,max( num_est-1, 1));
else 443 443 else
velo_avg(num_avgpkt_i) = CFO_v_estimate; 444 444 velo_avg(num_avgpkt_i) = CFO_v_estimate;
end 445 445 end
%average smoothing 446 446 %average smoothing
if num_avgpkt_i == 5 447 447 if num_avgpkt_i == 5
ant_est_velo(1,num_est) = mean(velo_avg) ; 448 448 ant_est_velo(1,num_est) = mean(velo_avg) ;
end 449 449 end
450 450
%derive h: function for theda estimation 451 451 %derive h: function for theda estimation
rx_cfo_corr_t = exp(-1i*2*pi*rx_cfo_est_final*[2176-1:2176+length(corr_target)-2]); 452 452 rx_cfo_corr_t = exp(-1i*2*pi*rx_cfo_est_final*[2176-1:2176+length(corr_target)-2]);
Gv_gin_idx = min( max( lts_peaks-512, 1) , length(raw_rx_dec) -1024 ); 453 453 Gv_gin_idx = min( max( lts_peaks-512, 1) , length(raw_rx_dec) -1024 );
rx_dec_cfo_corr = raw_rx_dec(Gv_gin_idx:Gv_gin_idx+1024-1).*rx_cfo_corr_t; 454 454 rx_dec_cfo_corr = raw_rx_dec(Gv_gin_idx:Gv_gin_idx+1024-1).*rx_cfo_corr_t;
455 455
%AoA estimation 456 456 %AoA estimation
v_est_walker = sqrt(100+v_obj^2); %m/s 457 457 v_est_walker = sqrt(100+v_obj^2); %m/s
lambda = v_c/center_freq; % wavelength 458 458 lambda = v_c/center_freq; % wavelength
delta = v_est_walker * sample_time; % spatial separation between successive antennas in the array 459 459 delta = v_est_walker * sample_time; % spatial separation between successive antennas in the array
% (twice the one-way separation to account for the round-trip time) 460 460 % (twice the one-way separation to account for the round-trip time)
Length_window = 15; 461 461 Length_window = 15;
462 462
yy = rx_dec_cfo_corr; 463 463 yy = rx_dec_cfo_corr;
xx = corr_target; 464 464 xx = corr_target;
Length_N = length(yy); 465 465 Length_N = length(yy);
hh = fft(yy)./fft(xx); 466 466 hh = fft(yy)./fft(xx);
A = zeros(19, Length_N); 467 467 A = zeros(19, Length_N);
A_max_ind = zeros(1, Length_N); 468 468 A_max_ind = zeros(1, Length_N);
A_temp = zeros(1, Length_N); 469 469 A_temp = zeros(1, Length_N);
theta_est = zeros(1, Length_N); 470 470 theta_est = zeros(1, Length_N);
471 471
for n = 1:(Length_N-Length_window) 472 472 for n = 1:(Length_N-Length_window)
473 473
for theta = (-pi/2):pi/18:(pi/2) 474 474 for theta = (-pi/2):pi/18:(pi/2)
475 475
theta_index = round((theta+pi/2)/(pi/18)+1); 476 476 theta_index = round((theta+pi/2)/(pi/18)+1);
477 477
for i0 = 1:Length_window 478 478 for i0 = 1:Length_window
A(theta_index,n) = A(theta_index,n) + hh(n+i0)*exp(1j*2*pi/lambda*i0*delta*sin(theta)); 479 479 A(theta_index,n) = A(theta_index,n) + hh(n+i0)*exp(1j*2*pi/lambda*i0*delta*sin(theta));
end 480 480 end
481 481
end 482 482 end
[ A_temp(n), A_max_ind(n)] = max(abs(A(:,n))); 483 483 [ A_temp(n), A_max_ind(n)] = max(abs(A(:,n)));
484 484
if A_max_ind(n)==1||A_max_ind(n)==19 485 485 if A_max_ind(n)==1||A_max_ind(n)==19