From a170391a9eff9d641bebb7f8723f809d9b9c3939 Mon Sep 17 00:00:00 2001 From: proto17 Date: Wed, 18 May 2022 23:30:40 -0400 Subject: [PATCH] Create test_integer_freq_offset.m --- .../test_integer_freq_offset.m | 215 ++++++++++++++++++ 1 file changed, 215 insertions(+) create mode 100644 matlab/updated_scripts/test_integer_freq_offset.m diff --git a/matlab/updated_scripts/test_integer_freq_offset.m b/matlab/updated_scripts/test_integer_freq_offset.m new file mode 100644 index 0000000..fab330d --- /dev/null +++ b/matlab/updated_scripts/test_integer_freq_offset.m @@ -0,0 +1,215 @@ +clear all; + +% This script is used to test time and frequency estimation / correction of +% a simulated DroneID burst. No equalization is done in this script. +% There is a simulated fractional time offset applied, but the script will +% fail if that time offset is anything other than 0 or 0.5. This is +% something that will likely have to be fixed by the equalizer or by +% upsampling to find the most correct fractional starting time offset + +%% Parameters +sample_rate = 15.36e6; +fft_size = get_fft_size(sample_rate); +carrier_spacing = round(sample_rate / fft_size); + +% Pick a worst case Parts Per Million (PPM) to support. 20 PPM is the +% value for the HackRF SDR's oscialltor. Then choose the worst case +% frequency that the SDR will operate at. In this case 5.9 GHz is the +% highest frequency the radio will need to tune to for DroneID. Use that +% to then figure out what the worst case frequency offset (in Hz) can be. +% That will be used as the search space for the Integer Frequency Offset +% (IFO) which is done later. +worst_case_ppm = 20; +worst_case_freq = 5.9e9; +max_allowed_freq_offset = worst_case_freq * worst_case_ppm / 1e6; +max_allowed_int_freq_offset = ceil(max_allowed_freq_offset / carrier_spacing); + +% Pick a frequency offset to test with +frequency_offset = max_allowed_freq_offset - 3.3e3; + +% Calculate the long and short cyclic prefix lengths +[long_cp_len, short_cp_len] = get_cyclic_prefix_lengths(sample_rate); + +% Below are the cyclic prefix lengths for each OFDM symbol +cyclic_prefixes = [ + 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 +]; + +% Calculate how many samples there are in one DroneID burst +total_burst_sample_count = (fft_size * length(cyclic_prefixes)) + sum(cyclic_prefixes); + +% Calculate how many guard carriers there are in one OFDM symbol +[left_guards, right_guards] = get_num_guard_carriers(sample_rate); + +% Number of constellation points in the underlying PSK +modulation_order = 4; + +% Number of occupied OFDM carriers per symbol +num_data_carriers = 600; + +%% OFDM Time Domain Creation +modulator = comm.OFDMModulator( ... + "FFTLength", fft_size, ... + "NumGuardBandCarriers", [left_guards; right_guards], ... + "CyclicPrefixLength", cyclic_prefixes, ... + "NumSymbols", length(cyclic_prefixes), ... + "InsertDCNull", true); + +% Modulate 9 random OFDM symbols worth of PSK samples +symbols = pskmod( ... + randi([0, modulation_order - 1], num_data_carriers, length(cyclic_prefixes)), ... + modulation_order); + +% Create the OFDM time domain +modulated_samples = modulator(symbols); + +%% Insert ZC Sequence at symbol 4 + +% Figure out where in the modulated samples that the ZC sequence should +% start (this includes the cyclic prefix) +symbol_4_offset = (fft_size * 3) + sum(cyclic_prefixes(1:3)); + +% Create the time domain ZC sequence for symbol 4 +zc_seq = create_zc(fft_size, 4); + +% Tack on the cyclic prefix to build the full symbol +full_zc_seq = [zc_seq(end-cyclic_prefixes(4)+1:end); zc_seq]; + +% Insert the full ZC sequence symbol in place of the modulated symbol 4 +modulated_samples(symbol_4_offset:symbol_4_offset + length(full_zc_seq) - 1) = full_zc_seq; + +% Convert to row vector +modulated_samples = reshape(modulated_samples, 1, []); + +%% Prepend/Append Zeros + +% These zeros are needed to simulate an actual received sample and the +% ambiguity of where the burst actually starts. The number of samples +% isn't critical, but should be something > about 30 samples +num_zeros = 50; + +zero_vec = zeros(1, num_zeros); +modulated_samples = [zero_vec, modulated_samples, zero_vec]; + +%% Add in a 0.5 sample time offset + +% This simulates the clocks of the transmitter and receiver not being +% perfectly synced up. Upsample by 4, then skip the first two samples, and +% downsample by 4 +modulated_samples = resample(modulated_samples, 4, 1); +modulated_samples = [modulated_samples(2:end), 0, 0]; +modulated_samples = resample(modulated_samples, 1, 4); + +%% Add Noise + +% This is not an exact SNR. Just add gaussian noise to the signal to +% simulate real world reception +modulated_samples = awgn(modulated_samples, 10, 'measured'); + +figure(1); +plot(10 * log10(abs(modulated_samples).^2)); +title('Time Domain with Noise'); +%% Apply Frequency Offset + +freq_offset_vector = exp(1j * 2 * pi * frequency_offset / sample_rate * [1:length(modulated_samples)]); +freq_offset_samples = modulated_samples .* freq_offset_vector; + +figure(2); +subplot(3, 1, 1); plot(10 * log10(abs(fftshift(fft(modulated_samples))).^2)); +title('Original Modulated Samples (FFT)') +subplot(3, 1, 2); plot(10 * log10(abs(fftshift(fft(freq_offset_vector))).^2)); +title('Frequency Offset Vector (FFT)') +subplot(3, 1, 3); plot(10 * log10(abs(fftshift(fft(freq_offset_samples))).^2)); +title('Frequency Offset Samples (FFT)') + +%% Find the ZC sequence + +% Find the start of the burst by looking at the cyclic prefixes. Make sure +% to look PAST `num_zeros` since in the real world it will not be known +% where the burst actually starts, so the search space will be rather large +sto_estimate = est_sto(freq_offset_samples, sample_rate, num_zeros * 1.5); + +% Calculate how many samples after the start of the burst that the ZC +% sequence should be located (this offset includes the ZC sequence cyclic +% prefix) +zc_seq_offset = (fft_size * 3) + long_cp_len + (short_cp_len * 3); + +% Calculate where in the sample vector the ZC sequence starts +auto_corr_index = sto_estimate + zc_seq_offset; + +fprintf("Sample offset error: %d\n", sto_estimate - num_zeros - 1); +%% Estimate coarse freq offset and adjust + +% Use the cyclic prefix of the ZC sequence to do coarse frequency offset +% estimation (the ZC sequence is the 4th OFDM symbol). +coarse_cp_len = cyclic_prefixes(4); +cp_len_window_size = coarse_cp_len; + +% The auto_corr_index points to the start of the ZC sequence, so back off +% by the length of the 4th cyclic prefix +coarse_offset_start = auto_corr_index - coarse_cp_len; + +% Extract out `fft_size + coarse_cp_len` samples +window = freq_offset_samples(coarse_offset_start:coarse_offset_start + fft_size + coarse_cp_len - 1); + +% Extract the first and last `coarse_cp_len` samples from the window +window_left = window(1:coarse_cp_len); +window_right = window(end - coarse_cp_len + 1:end); + +% Calculate the coarse frequency offset in radians and Hz by taking the dot +% product of the two windows +cfo_est_radians = angle(sum(window_left .* conj(window_right))) / fft_size; +cfo_est_hz = cfo_est_radians * sample_rate / (2 * pi); + +% Adjust for the frequency offset +cfo_est_adj_vector = exp(1j * cfo_est_radians * [1:length(freq_offset_samples)]); +freq_offset_samples = freq_offset_samples .* cfo_est_adj_vector; + +%% Integer frequency offset estimation +% The ZC sequence will be used to find the integer frequency offset. So, +% extract out just the samples for the ZC sequence (not including the +% cyclic prefix) +received_zc_sym = freq_offset_samples(auto_corr_index:auto_corr_index + fft_size - 1); + +figure(3); +plot(abs(received_zc_sym)); +title('ZC Sequence (Time Domain Magnitude)'); + +% Estimate the integer frequency offset +est_int_offset_hz = est_integer_freq_offset(received_zc_sym, sample_rate, max_allowed_int_freq_offset); + +% Adjust for the integer frequency offset +ifo_est_adj_vector = exp(1j * 2 * pi * est_int_offset_hz / sample_rate * [1:length(freq_offset_samples)]); +freq_offset_samples = freq_offset_samples .* ifo_est_adj_vector; + +% Print out how far off the combination of coarse and integer freq offset +% was from the true value +residual_freq_error_hz = est_int_offset_hz + cfo_est_hz + frequency_offset; +fprintf("Remaining Frequency Offset (Hz): %f\n", residual_freq_error_hz); + +%% Extract OFDM Data Carriers +demod = comm.OFDMDemodulator( ... + "FFTLength", fft_size, ... + "NumGuardBandCarriers", [left_guards; right_guards], ... + "CyclicPrefixLength", cyclic_prefixes, ... + "NumSymbols", length(cyclic_prefixes), ... + "RemoveDCCarrier", true); + +% Extract out just the samples in the burst and run those through the OFDM +% demodulator +burst_samples = freq_offset_samples(sto_estimate:sto_estimate + total_burst_sample_count - 1); +demod_symbols = demod(reshape(burst_samples, [], 1)); + +% Plot just the data carrying OFDM symbol constellations (the 4th symbol is +% the ZC sequence) +figure(4); +plot(demod_symbols(:,[1,2,3,5,6,7,8,9]), 'o'); +title('Demodulated Constellations')