mirror
/
dji_droneid
kopia lustrzana https://github.com/proto17/dji_droneid
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

*Lots* of changes 

These changes have been over several months and can't really be split apart meaningfully :\
main
David Protzman 2022-08-28 01:10:24 -04:00
rodzic fe761e9473
commit a731a12660
1 zmienionych plików z 88 dodań i 35 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

123
matlab/updated_scripts/process_file.m
Wyświetl plik

@ -45,6 +45,16 @@ filter_taps = fir1(filter_tap_count, signal_bandwidth/file_sample_rate); % Creat
%% Burst Extraction
[long_cp_len, short_cp_len] = get_cyclic_prefix_lengths(file_sample_rate);
cyclic_prefix_schedule = [
long_cp_len, ...
short_cp_len, ...
short_cp_len, ...
short_cp_len, ...
short_cp_len, ...
short_cp_len, ...
short_cp_len, ...
short_cp_len, ...
long_cp_len];
fft_size = get_fft_size(file_sample_rate);
% A correlation figure number of -1 will prevent plotting by the find_zc_indices_by_file function
@ -72,8 +82,7 @@ scrambler_x2_init = fliplr([0 0 1, 0 0 1 0, 0 0 1 1, 0 1 0 0, 0 1 0 1, 0 1 1 0,
% first symbol. Skipping the first symbol for those drones that have 9 OFDM symbols results in the new "first" symbol
% having a short cyclic prefix as well. So, since the burst extractor always assumes that there are 9 symbols, the
% first symbol is skipped for the purposes of coarse CFO. The second symbol is assumed to have a short cyclic prefix
coarse_cfo_symbol_sample_offset = fft_size + long_cp_len + 1;
coarse_cfo_symbol_cyclic_prefix = short_cp_len;
cfo_estimation_symbol_idx = 2;
%% Burst Processing
for burst_idx=1:size(bursts, 1)
@ -82,6 +91,7 @@ for burst_idx=1:size(bursts, 1)
if (enable_plots)
figure(43);
subplot(2, 1, 1);
plot(10 * log10(abs(burst).^2));
title('Time domain abs^2 10log10 (original)');
@ -147,36 +157,70 @@ for burst_idx=1:size(bursts, 1)
%% Apply low pass filter
burst = filter(filter_taps, 1, burst);
% Remove the extra samples at the front.
% TODO(15April2022) Honestly not sure why this needs to be 1.5, but it does...
offset = filter_tap_count * 1.5;
burst = burst(offset-1:end);
if (enable_plots)
figure(43);
subplot(2, 1, 2);
plot(10 * log10(abs(burst).^2));
title('Time domain abs^2 10log10 (filtered)')
end
%% Interpolate and find the true starting sample offset
interp_factor = 1;
burst = resample(burst, interp_factor, 1);
true_start_index = find_sto_cp(burst, file_sample_rate * interp_factor);
burst = resample(burst(true_start_index:end), 1, interp_factor);
% Plot cyclic prefixes overlayed with the replica from the end of the OFDM symbol
if (enable_plots)
offset = 1;
figure(7777);
for cp_idx=1:length(cyclic_prefix_schedule)
subplot(3, 3, cp_idx);
symbol = burst(offset:offset + cyclic_prefix_schedule(cp_idx) + fft_size - 1);
left = symbol(1:cyclic_prefix_schedule(cp_idx));
right = symbol(end - cyclic_prefix_schedule(cp_idx) + 1:end);
plot(abs(left));
hold on
plot(abs(right));
hold off;
title(['Cyclic Prefix Overlay ', mat2str(cp_idx)]);
offset = offset + length(symbol);
end
end
%% Coarse frequency offset adjustment using one of the OFDM symbols (see coarse_cfo_symbol_sample_offset definition)
% Get the cyclic prefix, and then the copy of the cyclic prefix that exists at the end of the OFDM symbol
cp = burst(...
coarse_cfo_symbol_sample_offset:...
coarse_cfo_symbol_sample_offset + coarse_cfo_symbol_cyclic_prefix - 1);
copy = burst(...
coarse_cfo_symbol_sample_offset + fft_size:...
coarse_cfo_symbol_sample_offset + fft_size + coarse_cfo_symbol_cyclic_prefix - 1);
% Get the expected starting index of the symbol to be used for CFO estimation
zc_start = long_cp_len + (fft_size * 3) + (short_cp_len * 3);
zc_start = zc_start + 6;
cfo_est_symbol = burst(zc_start - short_cp_len:zc_start + fft_size - 1);
% Get the cyclic prefix, and then the copy of the cyclic prefix that exists at the end of the OFDM symbol
cyclic_prefix = cfo_est_symbol(1:short_cp_len);
symbol_tail = cfo_est_symbol(end - short_cp_len + 1:end);
skip = 0;
cyclic_prefix = cyclic_prefix(skip+1:end-skip);
symbol_tail = symbol_tail(skip+1:end-skip);
% Calculate the frequency offset by taking the dot product of the two copies of the cyclic prefix and dividing out
% the number of samples in between each cyclic prefix sample (the FFT size)
offset_radians = angle(dot(cp, copy)) / fft_size;
offset_radians = angle(dot(cyclic_prefix, symbol_tail)) / fft_size;
offset_hz = offset_radians * file_sample_rate / (2 * pi);
if (enable_plots)
figure(999);
plot(abs(cp).^2);
plot(abs(cyclic_prefix).^2);
hold on;
plot(abs(copy).^2);
plot(abs(symbol_tail).^2, '*-', 'Color', 'red');
hold off;
title('Cyclic Prefix Overlay - CFO Estimate')
end
% Apply the inverse of the estimated frequency offset back to the signal
burst = burst .* exp(1j * -offset_radians * [1:length(burst)]);
burst = burst .* exp(1j * -offset_radians * [1:length(burst)]);
%% OFDM Symbol Processing
% Extract the individual OFDM symbols without the cyclic prefix for both time and frequency domains
@ -203,8 +247,10 @@ for burst_idx=1:size(bursts, 1)
figure(441);
subplot(2, 1, 1);
plot(abs(channel1).^2, '-');
title('ZC Sequence 1 Channel')
subplot(2, 1, 2);
plot(abs(channel2).^2, '-');
title('ZC Sequence 2 Channel')
end
% Only use the fisrt ZC sequence to do the initial equaliztion. Trying to use the average of both ends up with
@ -218,18 +264,12 @@ for burst_idx=1:size(bursts, 1)
% This is done for symbols 4 and 6 even though they contain ZC sequences. It's just to keep the logic clean
for idx=1:size(bits, 1)
% Equalize just the data carriers
data_carriers = freq_domain_symbols(idx,data_carrier_indices) .* channel;
data_carriers = freq_domain_symbols(idx,data_carrier_indices);
% Adjust for the walking phase offset that will be present if the first time domain sample wasn't sampled at
% just the right moment (fractional time offset). If there is any fractional time offset then in the freq
% domain there will be a phase offset that accumulates at each FFT bin. This causes a smearing that can be
% fixed by the channel estimation, but because there are no pilots the absolute phase is only correct for the
% OFDM symbols next to the symbol used for equalization. So, the absolute phase offset caused by the fractional
% time offset is adjusted by multiplying the phase offset by how far each OFDM symbol is from the one that was
% used to do equalization. Using symbol 5 because it's in the middle of the two ZC sequences, and so whatever
% phase offset was calculated between the two ZC's applies directly to OFDM symbol 5.
data_carriers = data_carriers .* exp(1j * (-channel_phase_adj * (idx - 5)));
if (enable_equalizer)
% Equalize just the data carriers
data_carriers = data_carriers .* channel;
end
% Demodulate/quantize the QPSK to bits
bits(idx,:) = quantize_qpsk(data_carriers);
@ -238,18 +278,31 @@ for burst_idx=1:size(bursts, 1)
figure(1);
subplot(3, 3, idx);
plot(data_carriers, 'o');
ylim([-1, 1]);
xlim([-1, 1]);
title(['Symbol ', mat2str(idx), ' IQ']);
figure(111);
subplot(3, 3, idx);
plot(10 * log10(abs(time_domain_symbols(idx,:)).^2), '-');
title(['Symbol ', mat2str(idx), ' Time Domain']);
figure(112);
subplot(3, 3, idx);
plot(10 * log10(abs(freq_domain_symbols(idx,:)).^2));
title(['Symbol ', mat2str(idx), ' Freq Domain']);
end
end
% Save the constellation plots to disk for debugging
% THIS CAN BE COMMENTED OUT IF NEEDED
saveas(gcf, sprintf('%s/images/ofdm_symbol_%d.png', this_script_path, burst_idx));
if (enable_plots)
% Save the constellation plots to disk for debugging
% THIS CAN BE COMMENTED OUT IF NEEDED
png_path = sprintf('%s/images/ofdm_symbol_%d.png', this_script_path, burst_idx);
try
saveas(gcf, png_path);
catch
error('Could not write out PNG file to "%s"', png_path);
end
end
% The remaining bits are descrambled using the same initial value, but more bits
second_scrambler = generate_scrambler_seq(7200, scrambler_x2_init);