From 95bc95e236f42d14c4a43e62d019723610f23024 Mon Sep 17 00:00:00 2001 From: David Protzman Date: Sun, 28 Aug 2022 01:07:46 -0400 Subject: [PATCH] Added integer frequency offset correction logic --- matlab/updated_scripts/process_file.m | 48 +++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/matlab/updated_scripts/process_file.m b/matlab/updated_scripts/process_file.m index cc30f26..319f479 100644 --- a/matlab/updated_scripts/process_file.m +++ b/matlab/updated_scripts/process_file.m @@ -100,6 +100,54 @@ for burst_idx=1:size(bursts, 1) grid on; end + %% Find Integer Frequency Offset + + % Exploiting the fact that during the first ZC sequence the DC carrier will be much lower in amplitude than the + % surrounding samples. Steps: + % 1. Extract just the time domain samples used in the first ZC sequence + % 2. Interpolate those time domain samples to increase the frequency resolution of the measurement + % 3. Get the power spectrum (abs squared of the FFT) + % 4. Look N elements around the center of the FFT for the lowest point (this is the center of the signal) + % 5. Calculate how far off from 0 Hz the lowest bin was, and frequency shift the upsampled signal by that value + % 6. Decimate the samples back to the original sample rate for further processing + + % Calculate the first sample index for the first ZC sequence (skipping the cyclic prefix) + offset = sum(cyclic_prefix_schedule(1:4)) + (fft_size * 3) + filter_tap_count; + + % Upsample (interpolate and filter) the ZC sequence samples + interp_rate = 10; + burst = resample(burst, interp_rate, 1); + + % Extract out just the samples for the first ZC sequence + zc_samples = burst((offset * interp_rate):(offset * interp_rate) + (fft_size * interp_rate) - 1); + + % Convert the time domain ZC sequence samples to the frequency domain + fft_bins = 10 * log10(abs(fftshift(fft(zc_samples))).^2); + + % Loop for the lowest bin in the middle of the frequency domain spectrum + bin_count = 15; % How far left and right to look for the lowest carrier + + % Set all of the FFT bins on the outside to infinity so they can't possibly be the minimum value + fft_bins(1:(fft_size * interp_rate / 2) - bin_count) = Inf; + fft_bins((fft_size * interp_rate / 2) + bin_count - 1:end) = Inf; + + % Find the index of the FFT bin with the lowest amplitude + [~, center_offset] = min(fft_bins); + + % Calculate the frequency needed to correct the integer offset, then conver that to radians + integer_offset = ((fft_size * interp_rate / 2) - center_offset + 1) * 15e3; + radians = 2 * pi * integer_offset / (file_sample_rate * interp_rate); + + % The initial radians measurement was taken well into the burst. Applying this as is would cause an absolute phase + % offset. The logic below is meant to prevent that problem. + radians_scalar = [0:length(burst) - 1];% - offset; + + % Apply a frequency adjustment + burst = burst .* exp(1j * radians * radians_scalar); + + % Downsample (filter and decimate) the burst samples + burst = resample(burst, 1, interp_rate); + %% Apply low pass filter burst = filter(filter_taps, 1, burst);