From a731a12660eb10d352d944287b49cddf6cee2b31 Mon Sep 17 00:00:00 2001 From: David Protzman Date: Sun, 28 Aug 2022 01:10:24 -0400 Subject: [PATCH] *Lots* of changes These changes have been over several months and can't really be split apart meaningfully :\ --- matlab/updated_scripts/process_file.m | 123 ++++++++++++++++++-------- 1 file changed, 88 insertions(+), 35 deletions(-) diff --git a/matlab/updated_scripts/process_file.m b/matlab/updated_scripts/process_file.m index 1f93b9b..d08a9b6 100644 --- a/matlab/updated_scripts/process_file.m +++ b/matlab/updated_scripts/process_file.m @@ -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);