Sayyed Dormiani-tabatabaei / Phase Corrected Channel Sounding - ECE 257B

Download as
Browse Code ยป

Commit 161bb9ffa3d2aecd2bdda585be82e44310b85cb6

Authored by Sayyed Dormiani-tabatabaei
0 parents
Exists in master

first commit

Showing 6 changed files with 463 additions and 0 deletions Side-by-side Diff

ESPRIT_ISA.m View file @ 161bb9f
... ... @@ -0,0 +1,29 @@
  1 +function [ ang_est, Es, En] = ESPRIT_ISA( Y, K, idx1, idx2)
  2 +%ESPRIT ALGORITHM ON DATA Y, ESTIMATING K ANGLES FROM ULA OF COHERENT SUBARRAYS OF
  3 +%TWO ELEMENTS. FIRST ELEMENTS OF EACH SUBARRAY AT INDEXES idx1, OTHER
  4 +%ELEMENTS IN idx2
  5 +% Y MEASUREMENTS
  6 +% K NUMBER OF DOAS
  7 +% idx1 LOCATION OF FIRST ELEMENT OF EACH SUBARRAY
  8 +% idx2 LOCATION OF SECOND ELEMENT OF EACH SUBARRAY
  9 +% ang_est ESTIMATE OF DOAS
  10 +% Es SIGNAL SUBSPACE
  11 +% En NOISE SUBSPACE
  12 +
  13 +% Ryy = cov(transpose(Y));
  14 +ang_est = zeros(K,1);
  15 +Ryy = Y*Y';
  16 +[E,D] = eig(Ryy);
  17 +[~,idx] = sort(diag(D),'descend');
  18 +for i = 1:K
  19 + Es = E(:,idx(i));
  20 + En = E(:,1:end-K);
  21 + Es1 = Es(idx1,:);
  22 + Es2 = Es(idx2,:);
  23 + Psi = pinv(Es1)*Es2;
  24 + tmp = (angle(eig(Psi)))/pi;
  25 + ang_est(i) = (asind((tmp)));
  26 +% ang_est(i) = tmp*90;
  27 +
  28 +end
  29 +end
RUSK_zc_sounder.mat View file @ 161bb9f

No preview for this file type

est_CFO.m View file @ 161bb9f
... ... @@ -0,0 +1,14 @@
  1 +function [cfo_est,cfo_ests] = est_CFO(IQ,N_STS,num_det)
  2 +%ESTIMATE CFO (CARRIER FREQUENCY OFFSET) FROM IQ DATA CONTAINING STS
  3 +%(SHORT TRAINING SEQUENCE) WHICH REPEATS num_det TIMES.
  4 +
  5 +IQ = reshape(IQ,1,length(IQ));
  6 +cfo_ests = zeros(num_det,1);
  7 +for j = 1:num_det
  8 + sts_ind1 = ( (j-1)*N_STS+1 ):( j*N_STS );
  9 + sts_ind2 = (sts_ind1(end)+1):sts_ind1(end)+N_STS;
  10 + cfo_ests(j) = IQ(sts_ind1)*IQ(sts_ind2)';
  11 +end
  12 +cfo_est = mean(angle(cfo_ests))/(pi*N_STS);
  13 +end
  14 +
gen_cal_sounder.m View file @ 161bb9f
... ... @@ -0,0 +1,45 @@
  1 +%% MIMO CALIBRATION SOUNDER
  2 +
  3 +clear all
  4 +clc
  5 +rng(1)
  6 +
  7 +%% PARAMETERS
  8 +N_TX = 8; %num transmitters
  9 +Tx_t = zeros(N_TX, 2^14);
  10 +N = 2^14;
  11 +
  12 +for i = 1:N_TX
  13 + F = zeros(1,N);
  14 + F(246+10*i) = exp(1i*2*pi*rand(1,1));
  15 + T = ifft(F,N);
  16 + Tx_t(i,:) = T./max(abs(T));
  17 +end
  18 +
  19 +Tx_t(2:7,:) = 0;
  20 +Tx_t = .8*[Tx_t(1,:)];
  21 +%% Write to binary files
  22 +for j = 1:1
  23 + filename1 = strcat('outside_real',num2str(j),'.bin');
  24 + filename2 = strcat('outside_imag',num2str(j),'.bin');
  25 + write_to_binary(Tx_t(j,:),filename1,filename2)
  26 +end
  27 +
  28 +% GBs = [1:N_Gs,(N/2)-N_Gs:(N/2)+N_Gs,N-N_Gs:N];%guard bands for sounder
  29 +% save('RUSK_zc_sounder.mat','sts_t','sts_f','zc','Tx_t','GBs','GI','F','T','Pilots','interval');
  30 +
  31 +%% Test processing
  32 +% rx = sum(Tx_t);
  33 +% det = abs(xcorr(rx,zc));
  34 +% det(det<200) = 0;
  35 +% det = det(length(rx)-1:end);
  36 +% [~,inds] = findpeaks(det);
  37 +% inds = sort(inds,'ascend');
  38 +% aoi = rx(inds(1)-1+length(zc):inds(1)-1+length(zc)+30*16);
  39 +%
  40 +% figure(2)
  41 +% plot(abs(xcorr(rx,zc)))
  42 +%
  43 +% figure(3)
  44 +% plot(abs(xcorr(aoi,sts_t(1,:))))
  45 +
prosses_test_8x12.m View file @ 161bb9f
... ... @@ -0,0 +1,365 @@
  1 +clear all
  2 +clc
  3 +rng(1)
  4 +%% Parameters
  5 +
  6 +load('RUSK_zc_sounder.mat')
  7 +
  8 +K = 8; %number DOAs
  9 +N_TX = 8; %num transmitters
  10 +N_RX = 12; %num receivers
  11 +N = 2^10; %number time samples in sounding sequence
  12 +L = 5; %snapshots per DOA estimation
  13 +DOA_subcarrier = 700; %subcarrier to perform DOA estimation
  14 +MAX_INDEX = 10e6; %stop index on rx data (to avoid processing GBs)
  15 +STRT_IDX = 1; %first data sample to observe
  16 +CFO_CORRECTION = 1; %toggle to correct for CFO
  17 +SAMPLE_CORRECTION = 1; %toggle to correct for sample start index
  18 +PHASE_CORRECTION = 1; %toggle to correct random phase offsets for each rx channel pair
  19 +PLOT_DET = 0; %plot toggle to see detections (use for debugging)
  20 +PLOT_STS = 0; %plot toggle to see STS detection (use for debugging)
  21 +PHASE_PLOT = 0; %plot toggle to see phase of each rx clock
  22 +CHANNEL_PLOT = 0; %plot toggle to observe the channel (use for debugging)
  23 +
  24 +%% Load received data
  25 +for j = 1:N_RX
  26 + filename = strcat('<path_to_data_folder>/RUSK_rx/RUSK_rx_ch',num2str(j),'.bin'); % edit this line! %get address of binary rx data files
  27 + fileID = fopen(filename); %open files
  28 + Q = fread(fileID,'float32','ieee-le'); %read files as IEEE little endian, float 32 format
  29 + a = Q(1:2:length(Q)); %inphase channel
  30 + b = Q(2:2:length(Q)); %quadrature channel
  31 + IQ_temp = [complex(double(a),double(b))]; %convert to complex numbers
  32 + IQ_temp = IQ_temp(STRT_IDX:end); %cut off beggining values
  33 + if j == 1
  34 + IQ = zeros(N_RX,min([length(IQ_temp),MAX_INDEX])); %if this is the first
  35 + end
  36 + if length(IQ_temp) >= size(IQ,2) %if it's too long shave some off
  37 + IQ_temp = IQ_temp(1:size(IQ,2));
  38 + end
  39 + if length(IQ_temp) < size(IQ,2) %if it's too short add some zeros
  40 + IQ_temp = [IQ_temp;zeros(size(IQ,2)-length(IQ_temp),1)];
  41 + end
  42 + IQ(j,:) = transpose(IQ_temp);
  43 + fclose(fileID);
  44 +end
  45 +clear IQ_temp a b filename fileID Q
  46 +
  47 +%% Detect starting frames
  48 +max_cap = round(size(IQ,2)/interval); %largest number of captures possible
  49 +inds = zeros(N_RX,max_cap); %allocate space for starting indices of detected signals
  50 +tx_ids = zeros(N_RX,max_cap); %allocate space for transmitter ID for each detected signal
  51 +cfos = zeros(N_RX,max_cap); %store CFO correction factor for each detection
  52 +
  53 +intervals = interval.*[0:3];
  54 +tmp = zeros(1,interval);
  55 +tmp(1:length(zc)) = zc;
  56 +zc_full = repmat(tmp,1,4);
  57 +
  58 +for k = 1:N_RX
  59 +
  60 + [detector,lags] = xcorr(IQ(k,:),zc); %detect signals by correlating rx data with ZC sequence
  61 + detector = abs(detector./max(abs(detector))); %normalize
  62 + d_thresh = 75*mean(detector); %arbitrary noise threshold
  63 + detector(detector<d_thresh) = 0; %shave off noise
  64 + detector(lags<=0) = 0; %fix errors from autocorr zero padding
  65 + detector = detector(size(IQ,2)-1:end-interval); %remove zero padded part of autocorr vector
  66 + [~,ind] = findpeaks(detector); %find detections
  67 + ind = sort(ind,'ascend'); %order detections by time of arrival
  68 + %% PLOT
  69 + if PLOT_DET
  70 + [detector2,lags] = xcorr(IQ(k,:),zc);
  71 + detector2 = abs(detector2./max(abs(detector2)));
  72 + detector2(lags<=0) = 0;
  73 + detector2 = detector2(size(IQ,2)-1:end-interval);
  74 + figure(1)
  75 + clf
  76 + subplot(2,1,1)
  77 + plot(abs(detector))
  78 + grid on
  79 + subplot(2,1,2)
  80 + plot(abs(detector2))
  81 + grid on
  82 + end
  83 +
  84 + %% Sort out multipath detections
  85 + tmp = zeros(1,length(ind)); %temporary array
  86 + region = 10; %number of samples on each side of detection to observe
  87 + for i = 1:length(ind) %for each packet detection
  88 + if ind(i)>0 %if we haven't already nulled this detection
  89 + val = ind(i); %store detection sample index
  90 + tmpind = abs(ind-val); %get relative displacement of all other detections
  91 + [~,candidates] = find(tmpind<region); %keep indicies of detections in the nearby region
  92 + [~,maxind] = max(detector(ind(candidates))); %find index of largest correlation, this is considered the true arrival
  93 + candidates(maxind) = []; %remove true arrival from candidates vector
  94 + ind(candidates) = 0; %set all remaining detections to zero
  95 + end
  96 + end
  97 + ind(ind==0) = []; %remove zeroed out detections
  98 + inds(k,1:length(ind)) = ind; %store in matrix to contain all detections over all rx channels
  99 +
  100 + %% Determine which Tx is talking for each packet detection
  101 + for i = 1:length(ind)
  102 + aoi = IQ(k,ind(i)+length(zc):ind(i) + length(zc)+45*length(sts_t)+1); %area of interest
  103 + det = sliding_corr(aoi,32); %autocorrelate
  104 + det = abs(det./max(abs(det))); %normalize
  105 + det(det<.7) = 0; %remove values below detection threshold
  106 + [~,tmp] = findpeaks(det); %find peaks
  107 + num_det = round(nnz(det)/length(sts_t)); %number of detections
  108 +
  109 + switch num_det %logic based on number of STS detected.
  110 + case {0,1,2,3,4} %Tx 1 sent 5 STS
  111 + tx_ids(k,i) = 1;
  112 + case {5,6,7,8,9} %Tx 2 sent 10 STS
  113 + tx_ids(k,i) = 2;
  114 + case {10,11,12,13,14} %Tx 3 sent 15 STS
  115 + tx_ids(k,i) = 3;
  116 + case {15,16,17,18,19} %Tx 4 sent 20 STS
  117 + tx_ids(k,i) = 4;
  118 + case {20,21,22,23,24} %Tx 5 sent 25 STS
  119 + tx_ids(k,i) = 5;
  120 + case {2526,27,28,29} %Tx 6 sent 30 STS
  121 + tx_ids(k,i) = 6;
  122 + case {30,31,32,33,34} %Tx 7 sent 35 STS
  123 + tx_ids(k,i) = 7;
  124 + case {35,36,37,38,39,40} %Tx 8 sent 40 STS
  125 + tx_ids(k,i) = 8;
  126 + end
  127 +
  128 + %% CFO estimation
  129 + if CFO_CORRECTION
  130 + [CFO,cfo_ests] = est_CFO(aoi,length(sts_t),num_det+1); %calculate CFO
  131 + cfos(k,i) = CFO; %record CFO
  132 + end
  133 + %% PLOT
  134 + if PLOT_STS
  135 + figure(2)
  136 + subplot(2,1,1)
  137 + plot(abs(det))
  138 + grid on
  139 + axis([-10,45*length(sts_t),0,1.15])
  140 + title('Autocorrelation at detection start index')
  141 + xlabel('Lag (samples)')
  142 + ylabel('Correlation magnitude')
  143 + subplot(2,1,2)
  144 + zplane(cfo_ests./mean(abs(cfo_ests)))
  145 + grid on
  146 + title('CFO estimates')
  147 + end
  148 + q=1;
  149 + end
  150 +end
  151 +
  152 +%% Go over detections, correct errors, insert detections below noise floor
  153 +
  154 +for j = 1:size(tx_ids,2) %for all detections
  155 + true_arrival = median(tx_ids(:,j)); %democratic vote on which arrival is correct
  156 + if true_arrival == 0 %if no arrival
  157 + last_idx = j-1; %this is the last informational index value in tx_ids, inds, and cfos
  158 + break %stop
  159 + end
  160 + for i = 1:N_RX %for each rx channel
  161 + est_arrival = tx_ids(i,j); %estimated tx channel
  162 + if est_arrival > true_arrival %when estimated tx channel > true tx channel
  163 + shift_right = est_arrival-true_arrival; %find how many spaces we'll need to shift right
  164 + if j == 1 %special logic for very first detection
  165 + tx_ids(i,:) = [(1:shift_right)+(true_arrival-1),tx_ids(i,1:end-shift_right)]; %shift tx_ids right by 'shift_right' number of indicies
  166 + inds(i,:) = [interval*(-1:-1*shift_right)+inds(i,1),inds(i,1:end-shift_right)]; %shift inds right by 'shift right' number of indicies
  167 + for k = 1:shift_right %deal with CFO estimation
  168 +% aoi = IQ(i, inds(i,k)+length(zc)+1 : inds(i,k)+length(zc)+45*length(sts_t)+1); %area of interest in rx samples
  169 +% [CFO,~] = est_CFO(aoi,length(sts_t),true_arrival*5-1); %estimate CFO
  170 + CFO = 0;
  171 + cfos(i,:) = [CFO,cfos(i,1:end-k)]; %shift current CFO values right by one
  172 + end
  173 + else %logic for general case (j ~= 1)
  174 + tx_ids(i,:) = [tx_ids(i,1:j-1),(1:shift_right)+(true_arrival-1),tx_ids(i,j:end-shift_right)]; %shift tx_ids right by 'shift_right' number of indicies
  175 + inds(i,:) = [inds(i,1:j-1), interval*(-1:-1*shift_right)+inds(i,j),inds(i,j:end-shift_right)]; %shift inds right by 'shift right' number of indicies
  176 + for k = 1:shift_right %deal with CFO estimation
  177 +% aoi = IQ(i, inds(i,j+k-1)+length(zc)+1 : inds(i,j+k-1)+length(zc)+45*length(sts_t)+1); %are of interest in rx samples
  178 +% [CFO,~] = est_CFO(aoi,length(sts_t),true_arrival*5-1); %estimate CFO
  179 + CFO = 0;
  180 + cfos(i,:) = [cfos(i,1:j+(k-1)-1),CFO,cfos(i,j:end-k)]; %shift current CFO values right by one
  181 + end
  182 + end
  183 + end
  184 + if true_arrival == N_TX && est_arrival < true_arrival %special logic for when true tx channel is channel N_TX, but estimated detection isn't
  185 + shift_right = est_arrival; %number of spaces to chift tx_ids, inds, and cfos right
  186 + tx_ids(i,:) = [tx_ids(i,1:j-1),(1:shift_right)+(true_arrival-1),tx_ids(i,j:end-shift_right)]; %shift tx_ids right by 'shift_right' number of indicies
  187 + inds(i,:) = [inds(i,1:j-1), interval*(-1:-1*shift_right)+inds(i,j),inds(i,j:end-shift_right)]; %shift inds right by 'shift_right' number of indicies
  188 + for k = 1:shift_right %deal with each individual CFO estimation
  189 +% aoi = IQ(i, inds(i,j+k-1)+length(zc)+1 : inds(i,j+k-1)+length(zc)+45*length(sts_t)+1); %are of interest in rx samples
  190 +% [CFO,~] = est_CFO(aoi,length(sts_t),true_arrival*5-1); %estimation CFO
  191 + CFO = 0;
  192 + cfos(i,:) = [cfos(i,1:j+(k-1)-1),CFO,cfos(i,j:end-k)]; %shift current CFO values right by one, insert newly calculated value
  193 + end
  194 + end
  195 + end
  196 +end
  197 +tx_ids = tx_ids(:,1:last_idx);
  198 +inds = inds(:,1:last_idx);
  199 +cfos = cfos(:,1:last_idx);
  200 +
  201 +%% Determine number of soundings per channel
  202 +num_soundings = zeros(N_RX,N_TX); %matrix to hold number of detections from specific Tx per rx channel
  203 +for i = 1:N_RX %for each Rx channel
  204 + for j = 1:N_TX %for each transmit channel
  205 + num_soundings(i,j) = sum(tx_ids(i,:) == j); %how many detections were from transmitter j
  206 + end
  207 +end
  208 +n_soundings = min(min(num_soundings)); %number of soundings heard
  209 +
  210 +
  211 +
  212 +
  213 +
  214 +%% ############################################################################
  215 +%% PROCESS DATA, FIND CHANNELS, CALCULATE DOA
  216 +%% ############################################################################
  217 +
  218 +%% For each transmitter, look at each receiver and calculate channel
  219 +H_ests = zeros(n_soundings,N,N_TX,N_RX); %allocate space for freq channel ests
  220 +h_ests = zeros(n_soundings,N,N_TX,N_RX); %allocate space for time channel ests
  221 +ang_ESPRIT = zeros(K,n_soundings-L); %allocate space for DOA estimates using ESPRIT
  222 +ang_APG = zeros(K,n_soundings-L); %allocate space for DOA estimates using APG
  223 +measurements = zeros(N_TX,N_RX,n_soundings); %allocate space for measurements
  224 +
  225 +for j = 1:N_TX
  226 + Rx_t = zeros(N_RX,N,n_soundings); %allocate space for time domain signals between each Rx and Tx j
  227 + for k = 1:N_RX %for each rx channel
  228 + [~, ind] = find(tx_ids(k,:) == j); %find start indices of soundings from Tx j on Rx k
  229 + for l = 1:length(ind) %for each sounding-
  230 + data_t = IQ(k,inds(k,ind(l))+length(zc)+45*length(sts_t)-1:inds(k,ind(l))+length(zc)+45*length(sts_t)+N-2); %get rx data
  231 +
  232 + %% Correct CFO
  233 + if CFO_CORRECTION
  234 + data_t = data_t.*exp(1i*2*pi*cfos(k,ind(l)).*[0:length(data_t)-1]);%correct CFO
  235 + end
  236 +
  237 + %% Process each data capture
  238 + data_f = fft(data_t,N); %fft of capture
  239 + H_est = data_f./transpose(Pilots); %channel estimate (freq)
  240 + H_est(GBs) = 0; %guard bands set to zero (otherwise divide by zero error)
  241 + H_ests(l,:,j,k) = H_est; %store channel estimate (freq)
  242 + h_est = ifft(H_est,N); %channel tap delays
  243 + h_ests(l,:,j,k) = h_est; %store channel tap delays
  244 + Rx_t(k,:,l) = H_est; %save raw rx data for DOA processing
  245 +
  246 + %% Check/Correct if start index was off by 1 or more samples
  247 + if SAMPLE_CORRECTION
  248 + [~,arrival_idx] = max(abs(h_est)); %peak detection of time domain channel response
  249 + if arrival_idx ~= 1 %if arrival wasn't on index 1 it must be corrected
  250 + [~,ind_error] = max(circshift(abs(h_est),N/2)); %find how far off the peak arrival was. Add circular shift to avoid wraparound errors
  251 + ind_error = ind_error-(N/2)-1; %correct for circular shift
  252 + data_t = IQ(k,inds(k,ind(l))+length(zc)+45*length(sts_t)-1+ind_error:inds(k,ind(l))+length(zc)+45*length(sts_t)+N-2+ind_error); %get data
  253 + if CFO_CORRECTION
  254 + data_t = data_t.*exp(1i*2*pi*cfos(k,ind(l)).*[0:length(data_t)-1]);%correct CFO
  255 + end
  256 + data_f = fft(data_t,N); %fft of capture
  257 + H_est = data_f./transpose(Pilots); %channel estimate (freq)
  258 + H_est(GBs) = 0; %guard bands set to zero (otherwise divide by zero error)
  259 + H_ests(l,:,j,k) = H_est; %store channel estimate (freq)
  260 + h_est = ifft(H_est,N); %channel tap delays
  261 + h_ests(l,:,j,k) = h_est; %store channel tap delays
  262 + end
  263 + Rx_t(k,:,l) = H_est; %save raw rx data for DOA processing
  264 + end
  265 +
  266 + %% PLOT
  267 + if CHANNEL_PLOT
  268 + figure(3)
  269 + subplot(2,1,1)
  270 + plot(abs(h_est))
  271 + axis([-10,N+10,0,1.25*max(abs(h_est))])
  272 + grid on
  273 + str = strcat('Channel Est, time, Rx ',num2str(k),', Tx ',num2str(j));
  274 + title(str)
  275 + xlabel('Time sample')
  276 + ylabel('Channel response magnitude')
  277 + pause(.01)
  278 + subplot(2,1,2)
  279 + plot(abs(H_est))
  280 + grid on
  281 + str2 = strcat('Channel Est, freq, Rx ',num2str(k),', Tx ',num2str(j));
  282 + title(str2)
  283 + xlabel('Subcarrier')
  284 + ylabel('Channel magnitude')
  285 + end
  286 + end
  287 + end
  288 +
  289 + for l = 1:n_soundings
  290 +
  291 + %% DOA estimation
  292 + data = squeeze(Rx_t(:,:,l));
  293 + measurements(j,:,l) = transpose(data(:,DOA_subcarrier));
  294 + if l > L
  295 + data = squeeze(measurements(j,:,l-(L-1):l));
  296 + [ang_est,~,~] = ESPRIT_ISA( data, 1,1:2:N_RX, 2:2:N_RX); %estimate DOA via ESPRIT
  297 + ang_ESPRIT(j,l-L) = ang_est; %store DOA
  298 + end
  299 + end
  300 +end
  301 +
  302 +if PHASE_CORRECTION
  303 + phase_offsets = zeros(N_TX,N_RX,N); %allocate space for rx channel phase correction
  304 + num_it = 100; %number of iterations
  305 + delta = 1; %number of snapshots to move per estimate
  306 + for j = 1:N_TX
  307 + Y = squeeze(measurements(j,:,:));
  308 + Ryy = Y*Y';
  309 + [U,S,V] = svd(Ryy);
  310 + Es = U(:,1);
  311 + A_phase = rand(N_RX,1);
  312 + for m = 1:num_it
  313 + A_current = exp(1i*(A_phase));
  314 + A_learned = Es*pinv(Es)*A_current;
  315 + argA = unwrap(angle(A_learned));
  316 + [U,S2,V] = svd(argA);
  317 + A_phase = U(:,1)*S2(1:1,1:1)*V(:,1:1)';
  318 + end
  319 + A_cal = exp(1i*(A_phase));
  320 + A_cal = A_cal.*repmat(conj(A_cal(1,:)),N_RX,1);
  321 + phase_offsets(j,:,i:i+delta-1) = repmat(A_cal(:,1),1,delta);
  322 + if PHASE_PLOT
  323 + a_vec = [];
  324 + for k = 1:N_RX
  325 + a_vec = [a_vec;ones(k,1)*A_cal(k,1)];
  326 + end
  327 + figure(6)
  328 + clf
  329 + zplane(a_vec)
  330 + axis(1.15.*[-1,1,-1,1])
  331 + str = strcat('capture ',num2str(i));
  332 + title(str)
  333 + pause(.001)
  334 + end
  335 + As = reshape(A_cal(:,1),2,N_RX/2);
  336 + Y= Y.*repmat(conj(A_cal(:,1)),1,size(Y,2));
  337 + for l = 1:n_soundings
  338 +
  339 + %% DOA estimation
  340 + if l > L
  341 + data = Y(:,l-(L-1):l);
  342 + [ang_est,~,~] = ESPRIT_ISA( data, 1,1:2:N_RX, 2:2:N_RX); %estimate DOA via ESPRIT
  343 + ang_ESPRIT(j,l-L) = ang_est; %store DOA
  344 + end
  345 + end
  346 + end
  347 +end
  348 +
  349 +
  350 +
  351 +%% Plot DOAs
  352 +WWII_badguys = [1,n_soundings-L,-25,25]; %axis values
  353 +figure(5)
  354 +clf
  355 +title('DOA estimates, ESPRIT, With Calibration')
  356 +hold on
  357 +plot(ang_ESPRIT')
  358 +grid on
  359 +xlabel('Capture number')
  360 +ylabel('DOA, \circ')
  361 +axis(WWII_badguys)
  362 +
  363 +
  364 +
  365 +
sliding_corr.m View file @ 161bb9f
... ... @@ -0,0 +1,10 @@
  1 +function [ auto_corr ] = sliding_corr(aoi,N)
  2 +%sliding autocorrelation
  3 +M = length(aoi);
  4 +aoi = reshape(aoi,1,M);
  5 +auto_corr = zeros(M-1,1);
  6 +for i = 1:length(aoi)-2*N-1
  7 + window = aoi(i:i+N-1);
  8 + auto_corr(i)= aoi(i+N:i+2*N-1)*window';
  9 +end
  10 +