Commit ca9f297f authored by Mark Ruzindana's avatar Mark Ruzindana
Browse files

Install script with no scalloping fix. Used for packet_gen.m right now.

parent 8efa733c
......@@ -7,12 +7,12 @@
#include <math.h>
#define BN_ELE 38 // Number of elements/antennas in the array
#define BN_BIN 20 // 25 // Number of frequency bins
#define BN_TIME 8000 // 4000 //40 // Number of decimated time samples
#define BN_BIN 25 // 20 // Number of frequency bins
#define BN_TIME 4000 // 8000 //40 // Number of decimated time samples
#define BN_BEAM 14 // Number of beams we are forming
#define BN_POL 4
#define BN_BEAM1 (BN_BEAM/2) // Number of beams we are forming
#define BN_TIME_STI 80 //40 // Number of decimated time samples per integrated beamformer output
#define BN_TIME_STI 40 //80 // Number of decimated time samples per integrated beamformer output
#define BN_STI (BN_TIME/BN_TIME_STI) // Number of short time integrations
#define BN_STI_BLOC 64
#define BN_ELE_BLOC 64
......
......@@ -7,9 +7,9 @@
#define NF 8 // Number of Fengines
#define NI 8 // Number of inputs per Fengine
#define NA (NF*NI) // Number of total antennas
#define NC 20 // 25 // Number of frequency channels
#define NC 25 // 20 // Number of frequency channels
#define NT 20 // Number of time samples per packet/mcnt
#define NM 400 // 200 // Number of packets/mcnts per block
#define NM 200 // 400 // Number of packets/mcnts per block
#define pow1 32768 // Next power of 2 >= Nm*Nt
#define nblocks2 (pow1/1024) // Block size for second kernel
#define nblocks1 (pow1/nblocks2) // Block size for first kernel
......
......@@ -60,21 +60,6 @@ __global__ void PFB_kernel(char2* pc2Data,
return;
}
// Discard channels and perform FFT shift (part of scalloping solution)
__global__ void Discard_Shift_kernel(float2* pf2FFTOut, float2* pf2DiscShift)
{
int pt = threadIdx.x; // N-point FFT index
int sb = blockIdx.x; // Number of elements x coarse channels (time series) index
int st = blockIdx.y; // Windows index (4000/32 = 125 windows)
int i = blockIdx.z; // Chunks of channels to recover
// Both pre-processor macros are defined in kernels.h
pf2DiscShift[fftshift_idx(pt,i,sb,st)].x = pf2FFTOut[recover_idx(pt,i,sb,st)].x;
pf2DiscShift[fftshift_idx(pt,i,sb,st)].y = pf2FFTOut[recover_idx(pt,i,sb,st)].y;
return;
}
// When PFB disabled just perform FFT.
__global__ void CopyDataForFFT(char2 *pc2Data, float2 *pf2FFTIn)
{
......
......@@ -4,17 +4,8 @@
#include <cuda.h>
#include <cufft.h>
#define N_POINTS_FFT 64
#define N_ELEMENTS 64
#define N_FINE_CHANS 5
// Discard channels and perform FFT shift
#define recover_idx(pt,i,sb,st) ((pt+(48*i)) + N_POINTS_FFT*(sb) + N_POINTS_FFT*N_ELEMENTS*N_FINE_CHANS*(st))
#define fftshift_idx(pt,i,sb,st) ((pt+(16*(1-i))) + (N_POINTS_FFT/2)*(sb) + (N_POINTS_FFT/2)*N_ELEMENTS*N_FINE_CHANS*(st))
// stuct of parameters for PFB. Values indicate default values.
//#define DEFAULT_PFB_PARAMS {4000, 32, 8, 25, 5, 64, 320, 0, (char*)"hanning\0", (char*)"float\0", (char*)"\0", 1};
#define DEFAULT_PFB_PARAMS {8000, 64, 8, 20, 5, 64, 320, 0, (char*)"hanning\0", (char*)"float\0", (char*)"\0", 1};
#define DEFAULT_PFB_PARAMS {4000, 32, 8, 25, 5, 64, 320, 0, (char*)"hanning\0", (char*)"float\0", (char*)"\0", 1};
// plot 1 mean to hide the plot of the filter before continuing.
typedef struct {
int samples;
......@@ -33,7 +24,6 @@ typedef struct {
__global__ void PFB_kernel(char2* pc2Data, float2* pf2FFTIn, float* pfPFBCoeff, params pfbParams);
__global__ void Discard_Shift_kernel(float2* pf2FFTOut, float2* pf2DiscShift);
__global__ void map(char* dataIn, char2* dataOut, int channelSelect, params pfbParams);
__global__ void CopyDataForFFT(char2* pc2Data, float2* pf2FFTIn);
__global__ void saveData(char2* dataIn, char2* dataOut);
......
......@@ -13,7 +13,6 @@ char2* g_pc2DataRead_d = NULL;
float2* g_pf2FFTIn_d = NULL;
float2* g_pf2FFTOut_d = NULL;
float2* g_pf2DiscShift_d = NULL;
float *g_pfPFBCoeff = NULL;
float *g_pfPFBCoeff_d = NULL;
......@@ -86,10 +85,6 @@ int runPFB(signed char* inputData_h, float* outputData_h, params pfbParams) {
while(!g_IsProcDone) {
//FFT
iRet = doFFT();
// New Code ///////////////////////////////////
Discard_Shift_kernel<<<g_dimGPFB, g_dimBPFB>>>(g_pf2FFTOut_d, g_pf2DiscShift_d);
/////////////////////////////////////////////////
if(iRet != EXIT_SUCCESS) {
(void) fprintf(stderr, "ERROR: FFT failed\n");
cleanUp();
......@@ -101,7 +96,6 @@ int runPFB(signed char* inputData_h, float* outputData_h, params pfbParams) {
// step input and output buffers.
g_pf2FFTIn_d += g_iNumSubBands * g_iNFFT;
g_pf2FFTOut_d += g_iNumSubBands * g_iNFFT;
g_pf2DiscShift_d += g_iNumSubBands * (g_iNFFT/2);
lProcData += g_iNumSubBands * g_iNFFT;
if(lProcData >= ltotData - NUM_TAPS*g_iNumSubBands*g_iNFFT){ // >= process 117 ffts leaving 256 time samples, > process 118 ffts leaving 224 time samples.
......@@ -120,14 +114,10 @@ int runPFB(signed char* inputData_h, float* outputData_h, params pfbParams) {
//wind back in/out ptrs - should put in another pointer as a process read ptr instead of updating the global ptr.
g_pf2FFTOut_d = g_pf2FFTOut_d - countFFT*g_iNumSubBands*g_iNFFT;
g_pf2FFTIn_d = g_pf2FFTIn_d -countFFT*g_iNumSubBands*g_iNFFT;
g_pf2DiscShift_d = g_pf2DiscShift_d - countFFT*g_iNumSubBands*(g_iNFFT/2);
//int outDataSize = countFFT * g_iNumSubBands * g_iNFFT;
// Modified variable outDataSize 1/2 of g_iNFFT due to the discard of half the channels //////
int outDataSize = countFFT * g_iNumSubBands * (g_iNFFT/2);
int outDataSize = countFFT * g_iNumSubBands * g_iNFFT;
//CUDASafeCallWithCleanUp(cudaMemcpy(outputData_h, fftOutPtr, outDataSize*sizeof(cufftComplex), cudaMemcpyDeviceToHost));
CUDASafeCallWithCleanUp(cudaMemcpy(outputData_h, g_pf2DiscShift_d, outDataSize*sizeof(cufftComplex), cudaMemcpyDeviceToHost));
////////////////////////////////////////////////////////////////////////////////////////////
CUDASafeCallWithCleanUp(cudaMemcpy(outputData_h, g_pf2FFTOut_d, outDataSize*sizeof(cufftComplex), cudaMemcpyDeviceToHost));
return iRet;
......@@ -312,20 +302,11 @@ int initPFB(int iCudaDevice, params pfbParams){
//int sizeDataBlock_in = g_iNumSubBands * g_iNFFT * sizeof(float2);
int sizeDataBlock_in = pfbParams.samples*g_iNumSubBands * sizeof(float2);
int sizeTotalDataBlock_out = pfbParams.samples*g_iNumSubBands * sizeof(float2); // output fft array same size as output data for convinence the full size is not used. In the pfb function the output data will be the fft counter times block amount in the fft.
// New variables ///////////////////////////////////
int g_iNwindows = 125; // with 8000 time samples, but with 4032 samples, Nwindows = 63
int sizeTotalDataBlock_disc = (g_iNFFT/2)*g_iNwindows*g_iNumSubBands * sizeof(float2); // Not sure about the value of this yet.
//////////////////////////////////////////////////
CUDASafeCallWithCleanUp(cudaMalloc((void **) &g_pf2FFTIn_d, sizeDataBlock_in));
CUDASafeCallWithCleanUp(cudaMalloc((void **) &g_pf2FFTOut_d, sizeTotalDataBlock_out)); // goal will be to update the output ptr each time it fires.
// New Code ////////////////////////////////////////////////////////
CUDASafeCallWithCleanUp(cudaMalloc((void **) &g_pf2DiscShift_d, sizeTotalDataBlock_disc)); // Discard and shift array
////////////////////////////////////////////////////////////////////
CUDASafeCallWithCleanUp(cudaMemset((void *) g_pf2FFTIn_d, 0, sizeDataBlock_in));
CUDASafeCallWithCleanUp(cudaMemset((void *) g_pf2FFTOut_d, 0, sizeTotalDataBlock_out));
// New Code ////////////////////////////////////////////////////////////
CUDASafeCallWithCleanUp(cudaMemset((void *) g_pf2DiscShift_d, 0, sizeTotalDataBlock_disc));
////////////////////////////////////////////////////////////////////////
// set kernel parameters
(void) fprintf(stdout, "\tSetting kernel parameters...\n");
......@@ -339,11 +320,8 @@ int initPFB(int iCudaDevice, params pfbParams){
g_dimGPFB.x = (g_iNumSubBands * g_iNFFT) / g_dimBPFB.x;
g_dimGCopy.x = (g_iNumSubBands * g_iNFFT) / g_dimBCopy.x;
g_dimGPFB.y = 125; // with 8000 time samples, but with 4032 samples, g_dimBPFB.y = 63
g_dimGCopy.y = 125; // same as g_dimGPFB.y
// g_dimBPFB.y = 63; // No. of windows given 64 point FFTs and 4000 time samples per block. Since 4000/64=62.5, I need to see whether 32 samples should be discarded.
// g_dimGPFB.y = 63; // No. of windows given 64 point FFTs and 4000 time samples per block
g_dimGPFB.y = 125;
g_dimGCopy.y = 125;
// map kernel params
mapGSize.x = pfbParams.samples;
......@@ -466,10 +444,6 @@ void cleanUp() {
(void) cudaFree(g_pf2FFTOut_d);
g_pf2FFTOut_d = NULL;
}
if (g_pf2DiscShift_d != NULL) {
(void) cudaFree(g_pf2DiscShift_d);
g_pf2DiscShift_d = NULL;
}
free(g_pfPFBCoeff);
(void) cudaFree(g_pfPFBCoeff_d);
......
......@@ -21,12 +21,12 @@
#define DEF_CUDA_DEVICE 0
#define DEF_SIZE_READ 262144 // data block size. should this be set dynamically once I get the data?
#define DEF_LEN_SPEC 64 // Transform size (previously 32)
#define DEF_LEN_SPEC 32 // Transform size
#define NUM_TAPS 8 // PFB Decimation factor
#define DEF_NUM_CHANNELS 20 // System spec for total number of channels (previously 25)
#define DEF_NUM_CHANNELS 25 // System spec for total number of channels
#define PFB_CHANNELS 5 // Number of coarse channels through PFB
#define DEF_NUM_ELEMENTS 64 // System spec for number of elements
#define SAMPLES 8000// Time samples. Used to be 4000 with 64 points
#define SAMPLES 4000// Time samples.
#define PFB_OUTPUT_BLOCK_SIZE (SAMPLES+3*32)*PFB_CHANNELS*DEF_NUM_ELEMENTS*2 // (3*DEF_LEN_SPEC is to add more samples on the end to make it look like 128 pfb windows had been processed for the pfb correlator)
......@@ -43,10 +43,6 @@
#define FILE_COEFF_DATATYPE "float"
#define FILE_COEFF_SUFFIX ".dat"
// // Discard channels and perform FFT shift
// #define recover_idx(pt,i,sb,st) ((pt+(47*i)) + DEF_LEN_SPEC*(sb) + DEF_LEN_SPEC*DEF_NUM_ELEMENTS*PFB_CHANNELS*(st))
// #define fftshift_idx(pt,i,sb,st) ((pt+(17*(1-i))) + DEF_LEN_SPEC*(sb) + DEF_LEN_SPEC*DEF_NUM_ELEMENTS*PFB_CHANNELS*(st))
#define USEC2SEC 1e-6
typedef unsigned char BYTE;
......
......@@ -13,13 +13,13 @@
#endif
#ifndef NFREQUENCY
//#define NFREQUENCY 25
#define NFREQUENCY 20
#define NFREQUENCY 25
//#define NFREQUENCY 20
#endif
#ifndef NTIME
//#define NTIME 4000
#define NTIME 8000
#define NTIME 4000
//#define NTIME 8000
#endif
#ifndef NTIME_PIPE
......
#define XGPU_VERSION 0.1+178@g61e94c6-dirty
#define XGPU_VERSION 0.1+179@g8efa733-dirty
......@@ -9,7 +9,7 @@
#include "config.h"
#define VERBOSE 0
#define SAVE 0
#define SAVE 1
#define DEBUG 0
// Total number of antennas (nominally 40)
......@@ -65,8 +65,8 @@
// Number of throwaway channels
// Changed variable for scalloping fix /////////////////////////////////
//#define N_CHAN_THROWAWAY 12
#define N_CHAN_THROWAWAY 112
#define N_CHAN_THROWAWAY 12
//#define N_CHAN_THROWAWAY 112
////////////////////////////////////////////////////////////////////////
// Total number of processed channels
......@@ -77,8 +77,8 @@
// Number of channels processed per XGPU instance?
// Changed variable for scalloping fix /////////////////////////////////
//#define N_CHAN_PER_X 25
#define N_CHAN_PER_X 20
#define N_CHAN_PER_X 25
//#define N_CHAN_PER_X 20
///////////////////////////////////////////////////////////////////////
#if N_CHAN_PER_X!=XGPU_NFREQUENCY
#warning "N_CHAN_PER_X must match frequency channels needed by xGPU"
......@@ -92,8 +92,8 @@
// Number of time samples processed per XGPU instance
// Changed variable for scalloping fix /////////////////////////////////
//#define N_TIME_PER_BLOCK 4000
#define N_TIME_PER_BLOCK 8000
#define N_TIME_PER_BLOCK 4000
//#define N_TIME_PER_BLOCK 8000
////////////////////////////////////////////////////////////////////////
#if N_TIME_PER_BLOCK!=XGPU_NTIME
#warning "N_TIME_PER_BLOCK must match the time samples needed by xGPU"
......
% Previously called interp_data.m
% Read correlations
clearvars;
%close all;
% White noise
% filename = '/lustre/pulsar/users/rprestag/FLAG/JUNK/JUNK/BF/2016_07_25_02:45:31A.fits';
% Correlated data
% filename = '/lustre/pulsar/users/rprestag/FLAG/JUNK/JUNK/BF/2016_07_25_02:37:42A.fits';
% Grid data
% filename = '/lustre/gbtdata/TGBT16A_508_01/BF/2016_07_26_23:45:30A.fits';
% info = fitsinfo(filename);
% data = fitsread(filename, 'binarytable', 1);
% tmp = data{3};
% R = zeros(2*length(tmp), 1);
% R(1:2:end) = real(tmp(1,:));
% R(2:2:end) = imag(tmp(1,:));
Nele = 40;
Nele_tot = 64;
Nbin = 25;
Nsamp = 4000;
Nbaselines_tot = (Nele_tot/2 + 1)*Nele_tot;
Nbaselines = (Nele + 1)*Nele/2;
Nblocks = (Nele_tot/2 + 1)*Nele_tot/4;
%big = figure();
%big2 = figure();
%small = figure();
tmp_idx = 1;
tmp = [1 6 11 17 22];
blk_rows = zeros(Nele_tot/2, Nele_tot/2);
for i = 1:Nele_tot/2
blk_rows(i,1:i) = (i-1)*i/2+1:(i-1)*i/2+i;
end
Rtot = zeros(Nele_tot, Nele_tot, Nbin);
% % PFB correlator path %%%%%%%%%%%%%%%%%%%%%%%%
% PATH = '/lustre/flag/';
% Coarse correlator path %%%%%%%%%%%%%%%%%%%%%%%%
PATH = '/lustre/flag/TGBT16A_508_01/TMP/BF/';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mcnt = [0];%, 200, 400, 600];
%for mcnt = 0:2:198
for k = 1:length(mcnt)
disp(['Processing mcnt=', num2str(mcnt(k))]);
FILE = fopen([PATH, sprintf('cor_mcnt_%d_B.out', mcnt(k))], 'r');
[R, count] = fscanf(FILE, '%g\n');
fclose(FILE);
for Nb = 1:Nbin
rb_real = R(2*Nbaselines_tot*(Nb - 1)+1:2:2*Nbaselines_tot*Nb);
rb_imag = R(2*Nbaselines_tot*(Nb - 1)+2:2:2*Nbaselines_tot*Nb);
rb = rb_real + 1j*rb_imag;
Rb = zeros(Nele_tot, Nele_tot);
for Nblk = 1:Nblocks
block_r = rb(4*(Nblk-1)+1:4*Nblk);
[row, col] = find(blk_rows == Nblk);
Rb(2*row - 1, 2*col - 1) = block_r(1);
if sum(diag(blk_rows) == Nblk) == 0
Rb(2*row - 1, 2*col) = block_r(2);
end
Rb(2*row , 2*col - 1) = block_r(3);
Rb(2*row , 2*col ) = block_r(4);
end
Rb = Rb + (Rb' - diag(diag(Rb'))); % Exploit symmetry
Rb = Rb./Nsamp;
Rtot(:,:,Nb) = Rtot(:,:,Nb) + Rb.*Nsamp;
% Commented out in recent file modified by mitch %%%%%%%%%%%%%
fig_mod = ceil(Nb/40);
fig_mod_plot = mod(Nb,40);
figure(fig_mod);
if fig_mod_plot == 0
fig_mod_plot = 40;
end
subplot(8,5,fig_mod_plot);
imagesc(abs(Rtot(1:Nele, 1:Nele, Nb)));
title(['Bin ', num2str(Nb)]);
drawnow;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
% % Added from more recent file to plot PFB output
% idx = 1:160;
% idx1 = reshape(idx, [5,32]);
% idx2 = idx1';
% stitch_idx = reshape(idx2, [160,1]);
%
% % stitch_idx = [];
% % for kk = 0:31
% % stitch_idx = [stitch_idx, kk:32:kk+128];
% % end
% % stitch_idx = stitch_idx + 1;
% % tmp2 = squeeze(Rtot(1,1,stitch_idx));
% %
% % mat_idx = [];
% % for jj = 0:4
% % mat_idx = [mat_idx, (jj:32:jj+128)'];
% % end
% % stitch_idx = [mat_idx; mat_idx+4; mat_idx+8; mat_idx+12; mat_idx+16; mat_idx+20; mat_idx+24] + 1;
% % keyboard;
% % tmp3 = squeeze(Rtot(1,1,stitch_idx));
%
% tmp2 = squeeze(Rtot(1,1,stitch_idx));
% figure(11);
% plot(0:length(tmp2)-1, 10*log10(abs(tmp2))); grid on;
% tmp = squeeze(Rtot(1,1,:));
% figure(10);
% plot(0:length(tmp)-1, 10*log10(abs(tmp))); grid on;
% Read correlations
clearvars;
%close all;
% White noise
% filename = '/lustre/pulsar/users/rprestag/FLAG/JUNK/JUNK/BF/2016_07_25_02:45:31A.fits';
% Correlated data
% filename = '/lustre/pulsar/users/rprestag/FLAG/JUNK/JUNK/BF/2016_07_25_02:37:42A.fits';
% Grid data
% filename = '/lustre/gbtdata/TGBT16A_508_01/BF/2016_07_26_23:45:30A.fits';
% info = fitsinfo(filename);
% data = fitsread(filename, 'binarytable', 1);
% tmp = data{3};
% R = zeros(2*length(tmp), 1);
% R(1:2:end) = real(tmp(1,:));
% R(2:2:end) = imag(tmp(1,:));
Nele = 40;
Nele_tot = 64;
Nbin = 160;
Nsamp = 125;
Nbaselines_tot = (Nele_tot/2 + 1)*Nele_tot;
Nbaselines = (Nele + 1)*Nele/2;
Nblocks = (Nele_tot/2 + 1)*Nele_tot/4;
%big = figure();
%big2 = figure();
%small = figure();
tmp_idx = 1;
tmp = [1 6 11 17 22];
blk_rows = zeros(Nele_tot/2, Nele_tot/2);
for i = 1:Nele_tot/2
blk_rows(i,1:i) = (i-1)*i/2+1:(i-1)*i/2+i;
end
Rtot = zeros(Nele_tot, Nele_tot, Nbin);
PATH = '/lustre/projects/flag/';
mcnt = [0]%, 200, 400, 600];
%for mcnt = 0:2:198
for k = 1:length(mcnt)
disp(['Processing mcnt=', num2str(mcnt(k))]);
FILE = fopen([PATH, sprintf('cor_mcnt_%d.out', mcnt(k))], 'r');
[R, count] = fscanf(FILE, '%g\n');
fclose(FILE);
for Nb = 1:Nbin
rb_real = R(2*Nbaselines_tot*(Nb - 1)+1:2:2*Nbaselines_tot*Nb);
rb_imag = R(2*Nbaselines_tot*(Nb - 1)+2:2:2*Nbaselines_tot*Nb);
rb = rb_real + 1j*rb_imag;
Rb = zeros(Nele_tot, Nele_tot);
for Nblk = 1:Nblocks
block_r = rb(4*(Nblk-1)+1:4*Nblk);
[row, col] = find(blk_rows == Nblk);
Rb(2*row - 1, 2*col - 1) = block_r(1);
if sum(diag(blk_rows) == Nblk) == 0
Rb(2*row - 1, 2*col) = block_r(2);
end
Rb(2*row , 2*col - 1) = block_r(3);
Rb(2*row , 2*col ) = block_r(4);
end
Rb = Rb + (Rb' - diag(diag(Rb'))); % Exploit symmetry
Rb = Rb./Nsamp;
Rtot(:,:,Nb) = Rtot(:,:,Nb) + Rb.*Nsamp;
% Commented out in recent file modified by mitch %%%%%%%%%%%%%
% fig_mod = ceil(Nb/40);
% fig_mod_plot = mod(Nb,40);
% figure(fig_mod);
% if fig_mod_plot == 0
% fig_mod_plot = 40;
% end
% subplot(8,5,fig_mod_plot);
% imagesc(abs(Rtot(1:Nele, 1:Nele, Nb)));
% title(['Bin ', num2str(Nb)]);
% drawnow;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
% Added from more recent file to plot PFB output
idx = 1:160;
idx1 = reshape(idx, [5,32]);
idx2 = idx1';
stitch_idx = reshape(idx2, [160,1]);
% stitch_idx = [];
% for kk = 0:31
% stitch_idx = [stitch_idx, kk:32:kk+128];
% end
% stitch_idx = stitch_idx + 1;
% tmp2 = squeeze(Rtot(1,1,stitch_idx));
%
% mat_idx = [];
% for jj = 0:4
% mat_idx = [mat_idx, (jj:32:jj+128)'];
% end
% stitch_idx = [mat_idx; mat_idx+4; mat_idx+8; mat_idx+12; mat_idx+16; mat_idx+20; mat_idx+24] + 1;
% keyboard;
% tmp3 = squeeze(Rtot(1,1,stitch_idx));
tmp2 = squeeze(Rtot(1,1,stitch_idx));
figure(11);
plot(0:length(tmp2)-1, 10*log10(abs(tmp2))); grid on;
tmp = squeeze(Rtot(1,1,:));
figure(10);
plot(0:length(tmp)-1, 10*log10(abs(tmp))); grid on;
% Read correlations
clearvars;
%close all;
% White noise
% filename = '/lustre/pulsar/users/rprestag/FLAG/JUNK/JUNK/BF/2016_07_25_02:45:31A.fits';
% Correlated data
% filename = '/lustre/pulsar/users/rprestag/FLAG/JUNK/JUNK/BF/2016_07_25_02:37:42A.fits';
% Grid data
% filename = '/lustre/gbtdata/TGBT16A_508_01/BF/2016_07_26_23:45:30A.fits';
% info = fitsinfo(filename);
% data = fitsread(filename, 'binarytable', 1);
% tmp = data{3};
% R = zeros(2*length(tmp), 1);
% R(1:2:end) = real(tmp(1,:));
% R(2:2:end) = imag(tmp(1,:));
Nele = 40;
Nele_tot = 64;
Nbin = 160;
Nsamp = 125;
Nbaselines_tot = (Nele_tot/2 + 1)*Nele_tot;
Nbaselines = (Nele + 1)*Nele/2;
Nblocks = (Nele_tot/2 + 1)*Nele_tot/4;
%big = figure();
%big2 = figure();
%small = figure();
tmp_idx = 1;
tmp = [1 6 11 17 22];
blk_rows = zeros(Nele_tot/2, Nele_tot/2);
for i = 1:Nele_tot/2
blk_rows(i,1:i) = (i-1)*i/2+1:(i-1)*i/2+i;
end
Rtot = zeros(Nele_tot, Nele_tot, Nbin);
PATH = '/lustre/projects/flag/';
mcnt = [0]%, 200, 400, 600];
%for mcnt = 0:2:198
for k = 1:length(mcnt)
disp(['Processing mcnt=', num2str(mcnt(k))]);
FILE = fopen([PATH, sprintf('cor_mcnt_%d.out', mcnt(k))], 'r');
[R, count] = fscanf(FILE, '%g\n');
fclose(FILE);
for Nb = 1:Nbin
rb_real = R(2*Nbaselines_tot*(Nb - 1)+1:2:2*Nbaselines_tot*Nb);
rb_imag = R(2*Nbaselines_tot*(Nb - 1)+2:2:2*Nbaselines_tot*Nb);
rb = rb_real + 1j*rb_imag;
Rb = zeros(Nele_tot, Nele_tot);
for Nblk = 1:Nblocks
block_r = rb(4*(Nblk-1)+1:4*Nblk);
[row, col] = find(blk_rows == Nblk);
Rb(2*row - 1, 2*col - 1) = block_r(1);
if sum(diag(blk_rows) == Nblk) == 0
Rb(2*row - 1, 2*col) = block_r(2);
end
Rb(2*row , 2*col - 1) = block_r(3);
Rb(2*row , 2*col ) = block_r(4);
end
Rb = Rb + (Rb' - diag(diag(Rb'))); % Exploit symmetry
Rb = Rb./Nsamp;
Rtot(:,:,Nb) = Rtot(:,:,Nb) + Rb.*Nsamp;
% Commented out in recent file modified by mitch %%%%%%%%%%%%%
% fig_mod = ceil(Nb/40);
% fig_mod_plot = mod(Nb,40);
% figure(fig_mod);
% if fig_mod_plot == 0
% fig_mod_plot = 40;
% end
% subplot(8,5,fig_mod_plot);
% imagesc(abs(Rtot(1:Nele, 1:Nele, Nb)));
% title(['Bin ', num2str(Nb)]);
% drawnow;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
% Added from more recent file to plot PFB output
idx = 1:160;
idx1 = reshape(idx, [5,32]);
idx2 = idx1';
stitch_idx = reshape(idx2, [160,1]);
% stitch_idx = [];
% for kk = 0:31
% stitch_idx = [stitch_idx, kk:32:kk+128];
% end
% stitch_idx = stitch_idx + 1;
% tmp2 = squeeze(Rtot(1,1,stitch_idx));
%
% mat_idx = [];
% for jj = 0:4
% mat_idx = [mat_idx, (jj:32:jj+128)'];
% end
% stitch_idx = [mat_idx; mat_idx+4; mat_idx+8; mat_idx+12; mat_idx+16; mat_idx+20; mat_idx+24] + 1;
% keyboard;
% tmp3 = squeeze(Rtot(1,1,stitch_idx));