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

Download as
Browse Code ยป

Commit e24d9d7ce6a4f5e60bfff39cecf3acf30ac76c89

Authored by Yujen Ku
1 parent 834ed0c586
Exists in master

add credit

Showing 1 changed file with 7 additions and 0 deletions Inline Diff

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