Small change in how the symbols are extracted

These changes have been over several months and can't really be split apart meaningfully :\

matlab/updated_scripts/extract_ofdm_symbol_samples.m

@@ -27,15 +27,16 @@ function [time_domain, freq_domain] = extract_ofdm_symbol_samples(samples, sampl
     sample_offset = 1;
     for idx=1:length(cp_lengths)
         % Skip the cyclic prefix
-        sample_offset = sample_offset + cp_lengths(idx);
+        symbol = samples(sample_offset:sample_offset + fft_size + cp_lengths(idx) - 1);
+        symbol = symbol(cp_lengths(idx) + 1:end);
 
         % Extract the time domain samples for this OFDM symbol
-        time_domain(idx,:) = samples(sample_offset:sample_offset+fft_size-1);
+        time_domain(idx,:) = symbol;
 
         % Convert the time domain samples into frequency domain
         freq_domain(idx,:) = fftshift(fft(time_domain(idx,:)));
         
-        sample_offset = sample_offset + fft_size;
+        sample_offset = sample_offset + fft_size + cp_lengths(idx);
     end
 end
 

matlab/updated_scripts/find_sto_cp.m

@@ -0,0 +1,65 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Find the start time offset based on the cyclic prefix
+%
+% This method is immune to the issues that plague the ZC sequences (frequency offset causes a time shift in the
+% correlation results)
+%
+% It's best to provide this function an upsampled copy of the burst to help fix any fractional time offset that might be
+% present
+%
+% @param samples Complex IQ samples that make up the full burst
+% @param sample_rate Sample rate (in Hz) of the provided samples
+% @return start_offset Sample index that the burst starts at (first sample of the first cyclic prefix)
+function [start_offset] = find_sto_cp(samples, sample_rate)
+    [long_cp_len, short_cp_len] = get_cyclic_prefix_lengths(sample_rate);
+    fft_size = get_fft_size(sample_rate);
+    cyclic_prefix_length_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];
+    num_ofdm_symbols = length(cyclic_prefix_length_schedule);
+
+    full_burst_len = sum(cyclic_prefix_length_schedule) + (fft_size * num_ofdm_symbols);
+    num_tests = length(samples) - full_burst_len;
+    scores_cp_sto = zeros(1, num_tests);
+
+    for idx=1:num_tests
+        offset = idx;
+        scores = zeros(1, num_ofdm_symbols);
+    
+        % Extract and correlate the samples that each cyclic prefix is expected
+        % to be at
+        for cp_idx=1:num_ofdm_symbols
+            cp_len = cyclic_prefix_length_schedule(cp_idx);
+    
+            % Extract the full OFDM symbol including cyclic prefix
+            window = samples(offset:offset + fft_size + cp_len - 1);
+    
+            % Extract the cyclic prefix and the final samples of the symbol
+            left = window(1:cp_len);
+            right = window(end - cp_len + 1:end);
+    
+            % Correlate the two windows
+            scores(cp_idx) = abs(xcorr(left, right, 0, 'normalized'));
+    
+            % Move the sample pointer forward by the full symbol size
+            offset = offset + cp_len + fft_size;
+        end
+    
+        % In the real DroneID the first OFDM symbol needs to be ignored since
+        % it isn't always present.  So, just average the correlation scores of
+        % all but the first element
+        scores_cp_sto(idx) = sum(scores(2:end)) / (length(scores) - 1);
+
+    end
+
+    % Find the index of the highest score
+    [~, start_offset] = max(scores_cp_sto);
+end
+

matlab/updated_scripts/process_file.m

@@ -27,15 +27,16 @@ if (~ isfile(turbo_decoder_path))
 end
 
 %% File Parameters
-enable_plots = 0; % Set to 0 to prevent the plots from popping up
-
-sample_type = 'single';
-file_path = '/opt/dji/collects/2437MHz_30.72MSPS.fc32';
-file_sample_rate = 30.72e6;
-file_freq_offset = 7.5e6; % This file was not recorded with the DroneID signal centered
-
+enable_plots = true;        % Set to false to prevent the plots from popping up
 correlation_threshold = 0.7; % The SNR is pretty good, so using a high correlation score (must be between 0.0 and 1.0)
 chunk_size = 10e6;           % Number of samples to process at a time
+enable_equalizer = true;     % Enable/disable the frequency domain equalizer
+
+%% Paramters that the user must change
+sample_type = 'single';
+file_path = 'YOUR_FILE_NAME_HERE';
+file_sample_rate = YOUR_SAMPLE_RATE_HERE;
+file_freq_offset = 0e6;
 
 %% Low Pass Filter Setup
 signal_bandwidth = 10e6; % The actual occupied bandwidth of the DroneID signal
@@ -44,11 +45,27 @@ 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
+correlation_fig_number = -1;
+if (enable_plots)
+    correlation_fig_number = 456;
+end
+
 % Making sure that the bursts that are extracted have enough padding for the low pass filter to start up and terminate
 bursts = extract_bursts_from_file(file_path, file_sample_rate, file_freq_offset, correlation_threshold, chunk_size,...
-    filter_tap_count, 'SampleType', sample_type, 'CorrelationFigNum', 456);
+    filter_tap_count, 'SampleType', sample_type, 'CorrelationFigNum', correlation_fig_number);
 
 assert(~isempty(bursts), "Did not find any bursts");
 
@@ -65,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)
@@ -75,42 +91,136 @@ 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)');
+        
+        % Plot the FFT, but average it with a single pole IIR filter to make it smoother
+        figure(1000);
+        fft_bins = 10 * log10(abs(fftshift(fft(burst))).^2);
+        running = fft_bins(1);
+        beta = 0.06;
+        for idx = 2:length(fft_bins)
+            running = (running * (1 - beta)) + (fft_bins(idx) * beta);
+            fft_bins(idx) = running;
+        end
+        x_axis = file_sample_rate / 2 * linspace(-1, 1, length(burst));
+        plot(x_axis, fft_bins);
+        title('Frequency Spectrum (averaged)');
+        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);
+    
+    % Apply a frequency adjustment
+    burst = burst .* exp(1j * radians * [0:length(burst) - 1]);
+    
+    % Downsample (filter and decimate) the burst samples
+    burst = resample(burst, 1, interp_rate);
+    
     %% 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
@@ -137,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
@@ -152,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);
@@ -172,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);