sweepRead.m

% Read in HackRF sweep data in binary format
%
% Gary Schafer, July 2026
%
% hackrf_sweep operates with a minimum 20 MHz steps.
%
% |********^********^********^********|********^********^********^********|
% |********^********^********^********|********^********^********^********|
% |********^********^********^********|********^********^********^********|
%
% where:
%
% |<->| = a step (20 MHz of spectrum)
% ^****^ = one block within a step (5 MHz of spectrum)
% Each line above = one sweep (consists of 1 or more steps)
%
% If the requested span is greater than 20 MHz, it will collect spectral
% data to cover spectrum at least equal to or greater than the requested
% span. For example, if the amount of spectrum request is 105 MHz, 
% hackrf_sweep will cover 6 steps since it would require at least 5 steps
% of 20 MHz; the remaining 5 MHz will be covered by an additional 20 MHz
% step, meaning a total of 6 steps.
% 
% Each 20 MHz step is divided into 5 MHz blocks. Based on how hackrf_sweep
% collects data, it will be recorded as 1-3-2-4. For example, a step covering
% 88 - 108 MHz (1 step) will record data in 4 blocks. The beginning frequency
% of each block will be, in order they are recorded, 88, 98, 93, and 103 MHz.

% The manner in which hackrf_sweep works is that:
% 1) It sets the sample rate to 20 MHz, but filters with a 15 MHz BW filter.
%    This can be seen in the verbose output of the hackrf_sweep command:
%    call hackrf_sample_rate_set(20.000 MHz)
%    call hackrf_baseband_filter_bandwidth_set(15.000 MHz)
% 2) It steps up 7.5 up from the low frequency.
% 3) It captures the energy within that 20 MHz span.
% 4) It calculates the FFT of the entire 20 MHz, but only uses half of it (-7.5 to -2.5 MHz &
%    2.5 to 7.5 MHz). NOTE: It windows the data with a Hann(ing) window.
% 5) It steps up 5 MHz and repeats step 4.
% 6) This gives you the data for one 20 MHz "step", but in the order of 1-3-2-4.
% 7) It steps up 20 MHz from the first frequency, and repeats steps 2 - 6.
%
% The output of the hackrf_sweep command, when set to binary format, is:
% 4 bytes: 32-bit unsigned integer with the number of bytes following for one 5-MHz-wide chunk.
%          This does not include this 4 byte value.
% 8 bytes: 64-bit unsigned integer with low frequency of 5 MHz step
% 8 bytes: 64-bit unsigned integer with high frequency of 5 MHz step
% X number of 4 bytes each: 32-bit floating point values of the spectral information (dBFS).
%                           "X" can be calculated from the 4 byte uint32 value by subtracting
%                           16 from it, then dividing by 4. For example, if the value of the
%                           initial 4 byte uint32 is 220, then the number of spectral values
%                           will be (220-16)/4 = 51.
%
% This worksheet will use three main variables.
% sweep = one sweep is the system tuning from the low frequency to the high frequency
% step = one step covers a 20 MHz chunk of spectrum. Each sweep will use an integer number of steps.
% block = each step consists of four (4) blocks of 5 MHz-wide spectral data.
%
% This worksheet assumes that each block is the same size (same number of spectral points), and each
% step consists of four (4) blocks.

clear all;

fid=fopen('FM-broadcast-sweep-50-900-5000BW-1000sweeps-20260716.bin','r');
file_struct=stat(fid);  % Read in data of the file
filesize=getfield(file_struct,'size');  % Read in the file size, in bytes.
block_bytes=fread(fid,1,'uint32')+4; % This is the number of bytes total in the first block
step_bytes=block_bytes*4;  % Calculate the total number of blocks in the step
total_blocks=filesize/(block_bytes);  % Calculate the total number of blocks in the file.
total_steps=total_blocks/4;  % Calculate the total number of 20 MHz steps in dataset.
pts=(block_bytes-20)/4; % This is the number of spectral points, subtracting the 4 bytes of the block
                      % size, the 16 bytes for the low and high frequency values, and dividing by 4
                      % to account for the 4 byte floating point values used for each spectral point.
f_low=fread(fid,1,'uint64'); % Read in first low frequency
f_high=fread(fid,1,'uint64'); % Read in the first high frequency

% Loop through each block to identify how many sweeps have occurred.

steps=1;
steps_per_sweep=1;

while(steps<=total_steps)

    fseek(fid,((steps*step_bytes)+4));
    f_temp=fread(fid,1,'uint32');
    if(f_temp==f_low)
        steps_per_sweep=steps;
        break
    endif
    steps=steps+1;
    
endwhile

f_max=f_low+(20e6*steps_per_sweep);
span=f_max-f_low;
sweep_bytes=steps_per_sweep*step_bytes;  % Calculate the number of bytes for each sweep.
sweep_total=total_steps/steps_per_sweep; % Calculate the total number of sweeps.

% Set basic variables

sweeps=0;
spec_ave=repmat([0],steps_per_sweep*pts*4,1);    % Setup array for creating averaged spectrum.
spec_max=repmat([-200],steps_per_sweep*pts*4,1); % Setup array for creating max spectrum.

% Start looping through each sweep

while(sweeps<sweep_total)

spec_data=[];  % Array for creating spectrum for each sweep.
steps=0;       % Track each step within a sweep.
blocks=1;      % Track each block within a step.
block_start=[0,2,1,3]; % This accounts for how each block is ordered within a step.

% Start looping through each step

while(steps<steps_per_sweep)

% Loop through each block within a step

    while(blocks<=4)
        current_point=sweeps*sweep_bytes+steps*step_bytes+(block_start(blocks))*block_bytes+20;
        fseek(fid,current_point);
        s=fread(fid,pts,'float');  % Read in 32-bit floats.
        spec_data=[spec_data;s];
        blocks=blocks+1;
    endwhile
steps=steps+1;
blocks=1;
endwhile
spec_ave=spec_ave+spec_data;      % Add up each sweep to average
spec_max=max(spec_max,spec_data); % Calculate the maximum on each sweep
sweeps=sweeps+1;
endwhile
fclose(fid);

spec_ave=spec_ave/sweep_total; % Normalize the average 

freq_pts=length(spec_ave);
span=f_max-f_low; % Calculate the span of the sweep.
f=((0:(freq_pts-1))*span/freq_pts)+f_low;
f=f';

% Display average and max traces

% Calculate the resolution bandwidth based on the Hann window.
% The Hann window has a NENBW of 1.5.
rbw=1.5*5e6/pts;
rbwText=['RBW: ' num2str(rbw) ' Hz'];

figure(1,'position',[100,100,1600,600]);
plot(f/1e6,spec_ave,'color','k','linewidth',2)
%line(f/1e6,spec_max)
xlim([f_low/1e6,f_max/1e6]);
grid on;
titleText=['Average Trace (' rbwText ')'];
title(titleText,'fontsize',24);
ylabel('Magnitude (dBFS)','fontsize',18);
xlabel('Frequency (MHz)','fontsize',18);

figure(2,'position',[100,100,1600,600]);
plot(f/1e6,spec_max,'color','k','linewidth',2)
xlim([f_low/1e6,f_max/1e6]);
grid on;
titleText=['Max Trace (' rbwText ')'];
title(titleText,'fontsize',24);
ylabel('Magnitude (dBFS)','fontsize',18);
xlabel('Frequency (MHz)','fontsize',18);