Added varargs for extract bursts and find zc seq funcs

matlab/updated_scripts/extract_bursts_from_file.m

@@ -4,6 +4,13 @@
 % Otherwise the start sample estimate from correlating for the ZC sequence will be off in time as well as the
 % correlation score being lower
 %
+% Function does accept the following varargs inputs:
+%
+% - SampleType: MATLAB numeric type that the samples in the provided file are stored as (ex: 'single', 'int16', etc)
+%               Defaults to 'single'
+% - CorrelationFigNum: Figure number to use for plotting the results of the cross correlation in find_zc_indices_by_file
+%                      Defaults to -1 which does not show the figure.  Valid values are -1, or > 0
+%
 % @param input_path File containing complex 32-bit floating point samples (interleaved I,Q,I,Q,...)
 % @param sample_rate Sample rate that the file was recorded at.  Must be an integer multiple of 15.36 MSPS (the minimum 
 %                    sample rate for the DroneID downlink)
@@ -16,15 +23,49 @@
 % @param padding How many additional samples before and after the burst to extract.  Must be >= 0
 % @return bursts A matrix where each row contains one burst
 function [bursts] = extract_bursts_from_file(input_path, sample_rate, frequency_offset, correlation_threshold,...
-    chunk_size, padding)
+    chunk_size, padding, varargin)
-
-    num_samples = get_sample_count_of_file(input_path);
     
-    lte_carrier_spacing = 15e3;                       % OFDM carrier spacing
+    assert(isstring(input_path) || ischar(input_path), "Input path must be a string or char array");
-    fft_size = sample_rate / lte_carrier_spacing;     % Number of samples per OFDM symbol (minus cyclic prefix)
+    assert(isnumeric(sample_rate), "Sample rate must be numeric");
-    long_cp_len = round(1/192000 * sample_rate);      % Number of samples in the long cyclic prefix
+    assert(sample_rate > 0, "Sample rate must be > 0")
-    short_cp_len = round(0.0000046875 * sample_rate); % Number of samples in the short cyclic prefix
+    assert(isnumeric(frequency_offset), "Frequency offset must be numeric");
+    assert(isnumeric(correlation_threshold), "Correlation threshold must be numeric");
+    assert(isnumeric(chunk_size), "Chunk size must be numeric");
+    assert(chunk_size > 0, "Chunk size must be > 0");
+    assert(isnumeric(padding), "Padding must be numeric");
+    assert(padding >= 0, "Padding must be >= 0");
+    assert(mod(length(varargin), 2) == 0, "Varargs length must be a multiple of 2");
+    
+    % Default the type of each I and Q value to 32-bit floating point
+    sample_type = 'single';
+    correlation_fig_num = -1;
+    
+    % Process the varargs inputs if they exist
+    for idx=1:2:length(varargin)
+        key = varargin{idx};
+        val = varargin{idx+1};
 
+        switch (key)
+            case 'SampleType'
+                sample_type = val;
+            case 'CorrelationFigNum'
+                correlation_fig_num = val;
+            otherwise
+                error('Invalid varargs key "%s"', key);
+        end
+    end
+
+    assert(isstring(sample_type) || ischar(sample_type), "SampleType must be a string or char array");
+    assert(isnumeric(correlation_fig_num), "CorrelationFigNum must be numeric");
+    assert(correlation_fig_num == -1 || correlation_fig_num > 0, "CorrelationFigNum must be -1 or > 0");
+    
+    % Get the number of complex IQ samples in the input file
+    num_samples = get_sample_count_of_file(input_path, sample_type);
+    
+    fft_size = get_fft_size(sample_rate);   % Number of samples per OFDM symbol (minus cyclic prefix)
+    [long_cp_len, short_cp_len] = get_cyclic_prefix_lengths(sample_rate);
+    
+    % Pre-calculate the frequency offset as a complex value
     freq_offset_constant = 1j * pi * 2 * (frequency_offset / sample_rate);
 
     % The first ZC sequence is the 4th symbol, and the `find_zc_indices_by_file` function will (assuming no major
@@ -33,7 +74,8 @@ function [bursts] = extract_bursts_from_file(input_path, sample_rate, frequency_
     zc_seq_offset = (fft_size * 4) + long_cp_len + (short_cp_len * 3);
     
     % Find all instances of the first ZC sequence
-    indices = find_zc_indices_by_file(input_path, sample_rate, frequency_offset, correlation_threshold, chunk_size);
+    indices = find_zc_indices_by_file(input_path, sample_rate, frequency_offset, correlation_threshold, chunk_size, ...
+        'SampleType', sample_type, 'CorrelationFigNum', correlation_fig_num);
 
     % In the DJI Mini 2 there are 9 OFDM symbols: 2 long cyclic prefixes, 7 short.  This isn't the case on all drones.
     % For some drones there are just 8 OFDM symbols.  It looks like those drones just don't send the first OFDM symbol

matlab/updated_scripts/find_zc_indices_by_file.m

@@ -1,42 +1,97 @@
-% This script only works for MATLAB right now.  Octave doesn't like something around the filtering section :(
+% Find all instances of the first ZC sequence in the provided file.  Does not have to read the entire file in at once
+%
+% Uses a cross correlator to search for the first ZC sequence in DroneID.
+%
+% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+% !!!  There is currently a bug in this script!  The correlation results are NOT normalized and might require the    !!!
+% !!!  the user to greatly increase or decrease the correlation threshold past the expected range of 0.0 - 1.0       !!!
+% !!!  In the code below there is a plotting function that is disabled by default.  Enabling it will plot the        !!!
+% !!!  correlation results so that you can better set the correlation treshold.                                      !!!
+% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+%
+% Varags:
+%   - SampleType: MATLAB numeric type that the samples in the provided file are stored as (ex: 'single', 'int16', etc)
+%                 Defaults to 'single'
+%   - CorrelationFigNum: Figure number to use when plotting the correlation results.  Defaults to not showing the
+%                        correlation results at all.  Must be > 0 to show the plot, or -1 to disable
+%
+% @param file_path Path to the complex IQ file to be processed
+% @param sample_rate Sample rate that the input file was recorded at
+% @param frequency_offset Amount of offset (in Hz, positive or negative) that the recording needs to be shifted by
+%                         before processing the samples.  This is useful when the recording was taken at an offset to
+%                         remove the DC spike
+% @param correlation_threshold Minimum correlation magnitude required to classify a score as containing a valid ZC
+%                              sequence.  Should be between 0.0 and 1.0.  Frequency offset (not accounted for by the
+%                              `frequency_offset` parameter) and low SNR can cause this value to need to be set lower
+%                              than normal.  A good starting value is 0.5
+% @param chunk_size Number of complex IQ samples to read into memory at a time.  Set this value as high as you can
+%                   without running out of memory.  The larger this buffer the faster the file will be processed
+% @param varargin (see above)
+% @return zc_indices A row vector containing the samples offsets where correlation scores were seen at or above the
+%                    specified threshold
+function [zc_indices] = find_zc_indices_by_file(file_path, sample_rate, frequency_offset, correlation_threshold, ...
+    chunk_size, varargin)
 
-%% Signal parameters
+    assert(isstring(file_path) || ischar(file_path), "Input file path must be a string or char array");
-% file_path = '/opt/dji/collects/2437MHz_30.72MSPS.fc32';
+    assert(isnumeric(sample_rate), "Sample rate must be numeric");
-% file_sample_rate = 30.72e6;     % Collected 2x oversampled
+    assert(sample_rate > 0, "Sample rate must be > 0")
-% rough_frequency_offset = 7.5e6; % The collected signal is 7.5 MHz off center
+    assert(isnumeric(frequency_offset), "Frequency offset must be numeric");
-% correlation_threshold = 0.7;    % Minimum correlation score to accept (0.0 - 1.0)
+    assert(isnumeric(correlation_threshold), "Correlation threshold must be numeric");
+    assert(isnumeric(chunk_size), "Chunk size must be numeric");
+    assert(chunk_size > 0, "Chunk size must be > 0");
+    assert(mod(length(varargin), 2) == 0, "Varargs length must be a multiple of 2");
+    
+    sample_type = 'single';
+    correlation_fig_num = -1;
 
+    for idx=1:2:length(varargin)
+        key = varargin{idx};
+        val = varargin{idx+1};
+
+        switch (key)
+            case 'SampleType'
+                sample_type = val;
+            case 'CorrelationFigNum'
+                correlation_fig_num = val;
+            otherwise
+                error('Invalid varargs key "%s"', key);
+        end
+    end
+
+    assert(ischar(sample_type) || isstring(sample_type), "SampleType must be a string or char array");
+    assert(isnumeric(correlation_fig_num), "CorrelationFigNum must be numeric");
+    assert(correlation_fig_num == -1 || correlation_fig_num > 0, "CorrelationFigNum must be -1, or > 0");
 
-function [zc_indices] = find_zc_indices_by_file(file_path, sample_rate, frequency_offset, correlation_threshold, chunk_size)
     %% LTE parameters
-    carrier_spacing = 15e3;
+    fft_size = get_fft_size(sample_rate);
     
     % Pre-calculate the frequency offset rotation
     freq_offset_constant = 1j * pi * 2 * (frequency_offset / sample_rate);
     
     % The output of the ZC sequence generator needs to be conjugated to be used in a filter.  Also create a variable that
     % can hold the filter state between chunks to prevent discontinuities
-    correlator_taps = conj(create_zc(sample_rate / carrier_spacing, 4));
+    correlator_taps = conj(create_zc(fft_size, 4));
     correlator_state = [];
     
+    % Figure out how many samples there are in the file
+    total_samples = get_sample_count_of_file(file_path, sample_type);
+    fprintf('There are %d samples in "%s"\n', total_samples, file_path);
+    
     % Open the IQ recording
     file_handle = fopen(file_path, 'r');
-    
+
-    % Figure out how many samples there are in the file
-    fseek(file_handle, 0, 'eof');
-    total_samples = ftell(file_handle) / 4 / 2; % 4 bytes per float, 2 floats per complex sample
-    fseek(file_handle, 0, 'bof');
-    fprintf('There are %d samples in "%s"\n', total_samples, file_path);
-    
     % Really large array to store the cross correlation results from *all* samples
     zc_scores = zeros(total_samples - length(correlator_taps), 1);
     
     sample_offset = 0;
     while (~ feof(file_handle))    
         %% Read in the next buffer
-        % The `fread` command will return interleaved real, imag values, so pack those into complex samples
+        % The `fread` command will return interleaved real, imag values, so pack those into complex samples making sure
-        floats = fread(file_handle, chunk_size * 2, 'float');
+        % that the resulting complex values are double precision (this is to prevent functions from complaining later)
-        samples = floats(1:2:end) + 1j * floats(2:2:end);
+        real_values = fread(file_handle, chunk_size * 2, sample_type);
+        samples = double(real_values(1:2:end) + 1j * real_values(2:2:end));
         
         %% Frequency shift the input
         % This is somewhat optional, but the correlation scores will go down fast if the offset is > 1 MHz
@@ -44,8 +99,7 @@ function [zc_indices] = find_zc_indices_by_file(file_path, sample_rate, frequenc
         samples = samples .* reshape(rotation_vector, [], 1);
         
         %% Correlate for the ZC sequence
-        % Use a FIR filter to search for the ZC sequence.  The resulting values will be normalized to between 0 and 1.0 if
+        % Use a FIR filter to search for the ZC sequence.  THIS WILL NOT NORMALIZE!!!
-        % the abs^2 is taken
         % TODO(9April2022): Would be nice to use the fftfilt function, but I don't know if it has issues with keeping state
         %                   as it doesn't have a state input like the filter command.  In testing the FFT filter is twice 
         %                   as fast
@@ -58,9 +112,11 @@ function [zc_indices] = find_zc_indices_by_file(file_path, sample_rate, frequenc
     % Get the floating normalized correlation results
     abs_scores = abs(zc_scores).^2;
     
-%     figure(1); 
+    if (correlation_fig_num > 0)
-%     plot(abs_scores);
+        figure(correlation_fig_num); 
-%     title('Correlation Scores (normalized)')
+        plot(abs_scores);
+        title('Correlation Scores (normalized)')
+    end
     
     % Find all places where the correlation result meets the specified threshold
     % This is going to find duplicates because there are very likely going to be two points right next to each other that

matlab/updated_scripts/process_file.m

@@ -47,7 +47,7 @@ fft_size = get_fft_size(file_sample_rate);
 
 % 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);
+    filter_tap_count, 'SampleType', sample_type, 'CorrelationFigNum', 456);
 
 assert(~isempty(bursts), "Did not find any bursts");