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/mission/status/ Breakthrough Listen data download: http://www.breakthroughinitiatives.org/OpenDataSearch 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');