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

Removed scripts that are no longer used

pull/2/head
David Protzman 2022-04-22 11:36:24 -04:00
rodzic 55a04c17b3
commit 9ab25d4d39
3 zmienionych plików z 0 dodań i 344 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

247
matlab/automatic_detector.m
Wyświetl plik

@ -1,247 +0,0 @@
clear all
%% Signal parameters
file_path = '/opt/dji/collects/2437MHz_30.72MSPS.fc32';
sample_rate = 30.72e6; % Collected 2x oversampled
critial_sample_rate = 15.36e6; % The actual sample rate for this signal
carrier_spacing = 15e3; % Per the LTE spec
signal_bandwidth = 10e6; % Per the LTE spec
rough_frequency_offset = 7.5e6; % The collected signal is 7.5 MHz off center
%% Read in the IQ samples for a single burst
% The burst is found using baudline and its cursor/delta time measurements
file_start_time = 0.802225;
file_start_sample_idx = uint64(file_start_time * sample_rate);
burst_duration_time = 0.000721;
burst_duration_samples = uint64(burst_duration_time * sample_rate);
samples = read_complex_floats(file_path, file_start_sample_idx, burst_duration_samples);
figure(1);
plot(10 * log10(abs(fftshift(fft(samples)))));
title('Original Samples')
%% Roughly center the signal of interest
% The signal might have been collected at an offset, so remove that offset
rotation_vector = exp(1j * 2* pi / sample_rate * (0:length(samples)-1) * rough_frequency_offset);
samples = samples .* rotation_vector;
figure(1000);
plot(10 * log10(abs(fftshift(fft(samples)))));
title('Original Samples - Freq Shifted')
orig_fft_size = sample_rate / carrier_spacing;
orig_long_cp_len = round(1/192000 * sample_rate);
orig_short_cp_len = round(0.0000046875 * sample_rate);
%% Low pass filter the original samples
filter_taps = fir1(200, signal_bandwidth/sample_rate);
samples = filter(filter_taps, 1, samples);
%% Search for the ZC sequence in symbol 4
% Create the ZC seqeunce
zc = create_zc(orig_fft_size, 4);
zc = reshape(zc, 1, []); % Reshape to match the samples vector
% Run a normalized cross correlation searching for the ZC sequence
% TODO(8Apr2022): Search through a smaller space (no need to look through everything)
scores = zeros(1, length(samples)-length(zc));
for idx=1:length(scores)
scores(idx) = normalized_xcorr(zc, samples(idx:idx + length(zc) - 1));
end
% Find the highest score value and index
[score, index] = max(abs(scores).^2);
if (score < 0.7)
warning("Correlation score for the first ZC sequence was bad. Might have errors later");
end
% The ZC sequence is the fourth symbol which means there are three FFT windows, one
% long cyclic prefix and two short cyclic prefixes before the ZC. But, to keep the
% correlation window as small as possible, the cyclic prefix of the ZC is ignored,
% so that adds an additional short cyclic prefix of offset
zc_seq_offset = (orig_fft_size * 3) + orig_long_cp_len + (orig_short_cp_len * 3);
% Calculate where the burst starts based on the correlation index
start_offset = round(index - zc_seq_offset);
angle(scores(index))
% samples = samples .* exp(1j * -angle(scores(index)) * (0:(length(samples)-1)));
samples = samples .* exp(1j * -angle(scores(index)));
% Trim all samples before the starting offset, and decimate at the same time.
% The decimation operation here is to get the signal back down to critical rate.
% Can get away with dropping every other sample because the signal is already filtered
samples = samples(start_offset:2:end);
%% Calculate critical rate parameters
fft_size = critial_sample_rate / carrier_spacing;
long_cp_len = round(0.0000052 * critial_sample_rate);
short_cp_len = round(0.00000469 * critial_sample_rate);
%% Plot symbol overlays
% The logic below is for debugging only. It draws OFDM symbol boundaries
% over a time domain view of the burst. It alternates between red and
% green transparent rectangles over the time domain as a sanity check that
% the time alignment is correct (starting at the right sample)
plot_symbol_boundaries(samples, critial_sample_rate, 4);
%% Coarse frequency adjustment
% My CFO detection isn't working, so rotate the constellation by hand (just eyeballing the constellation plots)
% offset_adj = -0.00005;
% samples = samples .* exp(1j * 2 * pi * offset_adj * (0:length(samples)-1));
%% Phase adjustment
% The phase angle adjustment made in the ZC sequence detection step got the QPSK to be lined up with points in
% [1,0], [0,1], [-1,0], and [0,-1] but we want it at [1,1], [-1,1], [-1,-1], and [1,-1] to make the demodulation
% decision simpler. That means that the constellation needs to be rotated by pi/2 (45 degrees).
% TODO(10Apr2022): Because of timing offsets there is a lot of "smearing" as the constellation rotates. Instead of
% fixing that right now, the constellation is rotated a little more to keep the points in distinct
% quadrants.
phase_offset = pi/2 - 0.5;
samples = samples * exp(1j * phase_offset);
%% Symbol extraction
symbols_time_domain = zeros(9, fft_size);
symbols_freq_domain = zeros(9, fft_size);
sample_offset = 1;
for symbol_idx=1:9
% OFDM symbols 1 and 9 have a long cyclic prefix
if (symbol_idx == 1 || symbol_idx == 9)
cyclic_prefix_len = long_cp_len;
else
cyclic_prefix_len = short_cp_len;
end
% Skip the cyclic prefix
sample_offset = sample_offset + cyclic_prefix_len;
% Extract out just the OFDM symbol ignoring the cyclic prefix
symbol = samples(sample_offset:sample_offset + fft_size - 1);
% Save off the original time domain and the FFT'd samples
symbols_time_domain(symbol_idx,:) = symbol;
symbols_freq_domain(symbol_idx,:) = fftshift(fft(symbol));
% Skip the OFDM symbol just extracted
sample_offset = sample_offset + fft_size;
end
% Plot the constellation diagrams for all extracted OFDM symbols
figure(5)
for idx=1:9
subplot(3, 3, idx);
plot(symbols_freq_domain(idx,:), 'o');
end
title('Constellations (Pre Freq Correction)')
%% Demodulation
% The spec says that the 15.36 MHz signal has 601 occupied carriers, but
% one of those is DC which doesn't appear to be used
% NOTE: I have seen someone saying that there is a specific phase present
% in some of the DC. Not sure if this is correct, so for now assuming that
% the DC carrier is not used
occupied_carriers = 600;
% Create a list of data carrier indices
data_carriers_left = (fft_size/2)-(occupied_carriers/2)+1:fft_size/2;
data_carriers_right = (fft_size/2)+2:fft_size/2+(occupied_carriers/2)+1;
% Below is logic to write the received ZC sequence out to a file. It is meant to be used with the
% brute_force_zc.m script
%
% % figure(1000); plot(abs(symbols_time_domain(4,:)).^2);
% file_handle = fopen('/opt/dji/collects/zc_symbol_4_15360KSPS_time_domain.fc32', 'w');
% zc_samples = symbols_freq_domain(4, [data_carriers_left, data_carriers_right]);
% reals = reshape([real(zc_samples); imag(zc_samples)], [], 1);
% fwrite(file_handle, reals, 'single');
% fclose(file_handle);
%
% % figure(1000); plot(abs(symbols_time_domain(6,:)).^2);
% file_handle = fopen('/opt/dji/collects/zc_symbol_6_15360KSPS_time_domain.fc32', 'w');
% zc_samples = symbols_freq_domain(6, [data_carriers_left, data_carriers_right]);
% reals = reshape([real(zc_samples); imag(zc_samples)], [], 1);
% fwrite(file_handle, reals, 'single');
% fclose(file_handle);
% Place to store the data carrier samples
data_carriers = zeros(7, occupied_carriers);
% Plot and extract the data carriers from the 7 data carrying OFDM symbols
figure(6);
output_idx = 1;
% Symbols 4 and 6 are ZC sequences and do not contain data
for idx=[1,2,3,5,7,8,9]
subplot(3, 3, idx);
% Get all of the freq domain samples for the current symbol, extract
% out just the data carriers, and plot
carriers = symbols_freq_domain(idx,:);
data_carriers(output_idx,:) = carriers([data_carriers_left, data_carriers_right]);
plot(squeeze(data_carriers(output_idx,:)), 'o');
output_idx = output_idx + 1;
end
% Plot the data carriers as one constellation (first one with lines, second
% with dots)
figure(7);
subplot(1, 2, 1);
plot(squeeze(data_carriers));
subplot(1, 2, 2);
plot(squeeze(data_carriers), 'o');
%% Demodulate the data carrier bits
% The first OFDM symbol is a constant value. It's a sync of sorts
% Initializing to 0x12345678 and reversing the bit order
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, 0 1 1 1, 1 0 0 0]);
% Generate enough random for one OFDM symbol (need two bits per data carrier since there are 2 bits per sample)
first_sym_scrambler_bit_count = length(data_carriers) * 2;
% Create enough scrambler to validate the first symbol as it will annihilate (XOR to all zeros) when the scrambler is
% applied.
first_sym_scrambler = generate_scrambler_seq(first_sym_scrambler_bit_count, scrambler_x2_init);
% Need to rotate the carriers by 180 degrees due to corrections that were made earlier
% TODO(13April2022): This needs to be removed as it's a correction that's required due to changes made earlier
data_carriers = data_carriers .* exp(1j * pi);
% Extract bits from the QPSK data carriers of the first OFDM symbol
first_data_carrier_bits = reshape(quantize_qpsk(data_carriers(1,:)), 1, []);
assert(isequal(bitxor(first_sym_scrambler, first_data_carrier_bits), zeros(size(first_sym_scrambler))),...
"Expected that the first data symbol would be all zeros after applying the scrambler...");
remaining_syms_scrambler = generate_scrambler_seq(6 * length(data_carriers) * 2, scrambler_x2_init);
bits = [];
for sym_idx = 2:7
bits = [bits; quantize_qpsk(data_carriers(sym_idx,:))];
end
descrambled_bits = bitxor(remaining_syms_scrambler, reshape(bits, 1, []));
handle = fopen("/tmp/bits", "wb");
fwrite(handle, descrambled_bits, 'int8');
fclose(handle);
zc2 = symbols_time_domain(6,:);
a = zc2(1:128);
scores = [];
for idx=1:(length(zc2)-length(a))
scores = [scores, normalized_xcorr(a, zc2(idx:idx+length(a)-1))];
end
figure(11111);
plot(abs(scores).^2);

59
matlab/brute_force_zc.m
Wyświetl plik

@ -1,59 +0,0 @@
% The input files for this script are generated by the `automatic_detector.m` script
clear all;
% Which symbol to select (must be 4 or 6 as those are the OFDM symbols that have zc sequences)
symbol_number = 4;
file_path = ['/opt/dji/collects/zc_symbol_', mat2str(symbol_number), '_15360KSPS_time_domain.fc32'];
file_handle = fopen(file_path, 'r');
reals = fread(file_handle, inf, 'single');
fclose(file_handle);
received_samples = reals(1:2:end) + 1j * reals(2:2:end);
figure(1); plot(10 * log10(abs(received_samples).^2));
root_sequnce_attempts = 1000;
scores = zeros(root_sequnce_attempts, 1);
for root=1:length(scores)
try
% Not all ZC roots are valid and the function will throw an error if the root is invalid
zc = zadoffChuSeq(root, 601);
% Remove the center sample as this is what's done on the transmitter.
zc = [zc(1:301); zc(303:end)];
catch
continue
end
% Correlate the
scores(root) = abs(xcorr(received_samples, fliplr(zc(1:end)), 0, 'normalized')).^2;
end
figure(2);
plot(scores);
title('ZC Normalized XCorr Results')
[value, index] = max(scores);
if (value < 0.7)
warning('Unable to find matching sequence');
end
zc = zadoffChuSeq(index, 601);
figure(3);
subplot(2, 1, 1);
spectrogram(zc);
title('Calculated ZC')
subplot(2, 1, 2);
spectrogram(received_samples);
title('Received ZC')
figure(4)
subplot(2, 1, 1);
plot(abs(fftshift(fft(zc))).^2);
title('Calculated ZC')
subplot(2, 1, 2);
plot(abs(fftshift(fft(received_samples))).^2);
title('Received ZC')

38
matlab/estimate_cp_freq_offset.m
Wyświetl plik

@ -1,38 +0,0 @@
function [freq_offset_est] = estimate_cp_freq_offset(samples, fft_size, short_cp_len, long_cp_len)
scores = zeros(9, 1);
sample_offset = 1;
figure(10);
for idx=1:9
if (idx == 1 || idx == 9)
cp_len = long_cp_len;
else
cp_len = short_cp_len;
end
symbol = samples(sample_offset:sample_offset + cp_len + fft_size - 1);
cyclic_prefix = symbol(1:cp_len);
end_of_symbol = symbol(end-cp_len+1:end);
% cyclic_prefix = samples(sample_offset+skip:sample_offset + cp_len - 1 - skip);
% end_of_symbol = samples(sample_offset+fft_size+skip:sample_offset+fft_size+cp_len - 1 - skip);
subplot(3, 3, idx);
plot(abs(cyclic_prefix).^2);
hold on;
plot(abs(end_of_symbol).^2);
hold off;
scores(idx) = angle(sum(cyclic_prefix .* conj(end_of_symbol))) / (2 * pi * cp_len);
% sum(angle(cyclic_prefix .* conj(end_of_symbol))) / length(fft_size) * 15.36e6
%scores(idx) = sum(angle(cyclic_prefix .* conj(end_of_symbol))) / length(cyclic_prefix);
sample_offset = sample_offset + fft_size + cp_len;
end
freq_offset_est = sum(scores) / 9;
scores * 15.36e6
%freq_offset_est = sum(scores) / length(scores) / fft_size;
end