Breakthrough Listen Raw Files Manipulation

Tutorial on Breakthrough Listen raw files handling for the detection of the Voyager 1 spacecraft.
Voyager missions:
https://www.jpl.nasa.gov/voyager/miss...
Breakthrough Listen data download:
http://www.breakthroughinitiatives.or...

MATLAB code:
fclose('all');
clear('all');
close('all');

nResol = 2^14; % resolution (in dObsbw/nObsnchan/nResol MHz)
NCHANOI = 7; % channel number

fname = 'blc07_guppi_57650_67573_Voyager1_0002.0005.raw'; % RAW file name

fid = fopen(fname,'r'); % open file

% read 1st header
header = fread(fid,80,'char')'; % instantiation of header
nHeaderLines = 1; % counter length of header
while ~strcmp(char(header(end,1:3)),'END') % until reaching keyword 'END'
header = [header;fread(fid,80,'char')'];
if strcmp(char(header(end,1:8)),'DIRECTIO')
nDirectio = str2num(char(header(end,10:end)));
end
if strcmp(char(header(end,1:7)),'OBSFREQ')
dObsfreq = str2num(char(header(end,10:end)));
end
if strcmp(char(header(end,1:8)),'BLOCSIZE')
nBlocsize = str2num(char(header(end,10:end)));
end
if strcmp(char(header(end,1:5)),'NBITS')
nBits = str2num(char(header(end,10:end)));
end
if strcmp(char(header(end,1:8)),'OBSNCHAN')
nObsnchan = str2num(char(header(end,10:end)));
end
if strcmp(char(header(end,1:5)),'OBSBW')
dObsbw = str2num(char(header(end,10:end)));
end
nHeaderLines = nHeaderLines+1;
end
nChanSize = nBlocsize/nObsnchan; % size of one channel in # of samples
nPadd = 0;
if nDirectio
nPadd = (floor(80*nHeaderLines/512)+1)*512 - 80*nHeaderLines; % directIO padding
end

fLow = dObsfreq - dObsbw/2; % lowest frequency of the band
fHigh = dObsfreq + dObsbw/2; % highest frequency of the band
dChanBW = dObsbw/nObsnchan; % bandwidth of single coarse channel

s = dir(fname); % size of whole file in bytes
NumBlocs = round(s.bytes/(nBlocsize+nPadd+80*nHeaderLines)); % number of data blocs

spec1 = zeros(nResol,NumBlocs);spec2 = zeros(nResol,NumBlocs); % instantiation of dynamic spectra
for nBLOC = 1:NumBlocs

fseek(fid, (nBLOC-1)*(nHeaderLines*80+nPadd+nBlocsize)+...
nHeaderLines*80+nPadd+(NCHANOI-1)*nChanSize, 'bof'); % go to right spot in file
data = fread(fid,nChanSize,'int8'); % read samples
data = reshape(data,2,length(data)/2); % separate real and imaginary
data = data(1,:) + 1i*data(2,:);
data = reshape(data,2,length(data)/2); % separate polarizations

% polarization #1
dataFT1 = reshape(data(1,1:floor(size(data,2)/nResol)*nResol),...
nResol,floor(size(data,2)/nResol));
dataFT1 = abs(fft(dataFT1,[],1)).^2; % averaged periodogram

% polarization #2
dataFT2 = reshape(data(2,1:floor(size(data,2)/nResol)*nResol),...
nResol,floor(size(data,2)/nResol));
dataFT2 = abs(fft(dataFT2,[],1)).^2;

spec1(:,nBLOC) = fftshift(mean(dataFT1,2));
spec2(:,nBLOC) = fftshift(mean(dataFT2,2));

waitbar(nBLOC/NumBlocs);
end

% lowest and highest frequencies of processed channel
fLowChan = fLow + (NCHANOI-1)*dChanBW;
fHighChan = fLowChan + dChanBW;

figure('Units','normalized','outerposition',[0,.4,1,.5]);
subplot(1,2,1);
imagesc(linspace(fLowChan,fHighChan,nResol),1:NumBlocs,10*log10(abs(spec1).^2)');
xlabel('frequency [MHz]');ylabel('block number');
title(['POLARIZATION #1 - integration:' num2str(NumBlocs*nChanSize/4/abs(dChanBW)/(10^6)) 's']);
subplot(1,2,2);
imagesc(linspace(fLowChan,fHighChan,nResol),1:NumBlocs,10*log10(abs(spec2).^2)');
xlabel('frequency [MHz]');ylabel('block number');
title(['POLARIZATION #2 - integration:' num2str(NumBlocs*nChanSize/4/abs(dChanBW)/(10^6)) 's']);

fclose('all');