Commit 6118582a authored by Mark Ruzindana's avatar Mark Ruzindana
Browse files

Updated MATLAB code

parent 3d58febf
% Simple script that saves the weights to a RTBF-compatible set of
% files
% Reshape weights
N_beam = 7;
N_ele = 64;
N_bin_total = 500;
N_bin = 25;
N_pol = 2;
% w_padded = zeros(N_ele, N_beam, N_pol, N_bin_total);
% w_padded(X_idx,:,1,:) = wX;
% w_padded(Y_idx,:,2,:) = wY;
w_padded = ones(N_ele, N_bin_total, N_beam, N_pol);
% Save data into weight file formatted for RTBF code
banks = {'A', 'B', 'C', 'D',...
'E', 'F', 'G', 'H',...
'I', 'J', 'K', 'L',...
'M', 'N', 'O', 'P',...
'Q', 'R', 'S', 'T'};
interleaved_w = zeros(2*N_ele*N_bin*N_beam*N_pol,1);
chan_idx = [1:5, 101:105, 201:205, 301:305, 401:405];
for b = 1:length(banks)
% Get bank name
bank_name = banks{b};
% Extract channels for bank
% w1 = w_padded(:,:,:,chan_idx+5*(b-1));
w1 = w_padded(:,chan_idx+5*(b-1),:,:);
%%% Just for testing - Used .mat file to verify that the .bin file was
%%% the same
% wp_filename = strrep(filename, '.bin', sprintf('_%s.mat',bank_name));
% save(wp_filename, 'w1');
% Reshape for file format
w2 = reshape(w1, N_ele*N_bin, N_beam*N_pol);
w_real = real(w2(:));
w_imag = imag(w2(:));
interleaved_w(1:2:end) = w_real(:);
interleaved_w(2:2:end) = w_imag(:);
% Get filename
filename = '/lustre/projects/flag/weight_files/dummy_w.bin';
weight_file = strrep(filename, '.', sprintf('_%s.', bank_name));
% Create metadata for weight file
offsets_el = 0.7;
offsets_az = 0.7;
offsets = [offsets_el; offsets_az];
offsets = offsets(:);
cal_filename = 'test';
to_skip1 = 64 - length(cal_filename);
algorithm_name = 'Max Signal-to-Noise Ratio';
to_skip2 = 64 - length(algorithm_name);
xid = b-1;
% Write to binary file
WID = fopen(weight_file,'wb');
if WID == -1
error('Author:Function:OpenFile', 'Cannot open file: %s', weight_file);
end
% Write payload
fwrite(WID, single(interleaved_w), 'single');
% Write metadata
fwrite(WID,single(offsets),'float');
fwrite(WID,cal_filename, 'char*1');
if to_skip1 > 0
fwrite(WID, char(zeros(1,to_skip1)));
end
fwrite(WID,algorithm_name, 'char*1');
if to_skip2 > 0
fwrite(WID, char(zeros(1,to_skip2)));
end
fwrite(WID, uint64(xid), 'uint64');
fclose(WID);
% fprintf('Saved to %s\n', weight_file);
end
\ No newline at end of file
......@@ -37,8 +37,9 @@ sigma2 = kb*Tsys*BW; % Noise power per channel
% 6 -> Send complex sinusodial data
% 7 -> Send ULA data
% 8 -> Send exponentially correlated noise.
% 9 -> Send pulsar data
% else -> Send all zeros
data_flag = 8;
data_flag = 9;
% Sinusoid parameters (only used if data_flag = 2)
% It should be noted that the phase of the sinusoid will not change between
......@@ -151,13 +152,43 @@ c_min = -4;
CEN_real = int8(((real(CEN) - c_min)/(c_max - c_min) - 0.5) * 256);
CEN_imag = int8(((imag(CEN) - c_min)/(c_max - c_min) - 0.5) * 256);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Case 9 - Simulated pulsar data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Increase the range of tau when dispersion measure causes m_D to exceed
% time samples.
D = 10; % Dispersion measure
freq = (0:499)*(303e3) + 1300e6; % All frequencies
fo = freq(floor(length(freq)/2)); % Center frequency
m_D = 4.1488e3*((fo^-2)-(freq.^-2))*D; % Frequency dependent timing offset
Ntime = 4000;
tau = -2.8e-14:((2.5e-14)+(2.8e-14))/(Ntime-1):2.5e-14; % Range of timing offsets
pulseData = zeros(Ninputs, Nbins, Ntime);
% Noisy environment
for ii = 1:size(pulseData,3)
for jj = 1:size(pulseData,2)
pulseData(:,jj,ii) = 0.1*(randn(1) + 1j*rand(1));
end
end
% Pulsar
pulse = 1;
for cyc = [0, -1000, 1000]
for m = 1:Ninputs
for k = 1:Nbins
tmp = abs(m_D(k)-tau);
[~,idx] = min(tmp);
phi = m*2*pi*freq(k)*tau(idx+cyc);
pulseData(m,k,idx+cyc) = pulse*exp(1j*phi) + 0.1*(randn(1) + 1j*rand(1));
end
end
end
% Create UDP sockets - 1 IP address per Xengine (xid)
for xid = 1:Nxengines
remoteHost = ['10.10.1.', num2str(xid)];
if xid == 1
remoteHost = '10.17.16.208';
remoteHost = '10.17.16.208'; % It was 208 before
end
%if xid == 14
% remoteHost = '10.10.1.1';
......@@ -172,7 +203,7 @@ end
mcnt = 0; % Each mcnt represents 20 packets across all F-engines in the
% same time frame
for mcnt = [0:201, 400:50:20000] %while mcnt <= 10000
for mcnt = [0:401]% [0:201, 400:50:20000] %while mcnt <= 10000
disp(['Sending mcnt = ', num2str(mcnt)]);
for xid = 1:1 % Set to a single X-engine for single HPC testing (Richard B.)
for fid = 1:Nfengines
......@@ -212,7 +243,7 @@ for mcnt = [0:201, 400:50:20000] %while mcnt <= 10000
w_idx = w_idx + 1;
% Generate signal to send
data = zeros(Nin_per_f, 2, Nbin_per_x, Ntime_per_packet);
data = zeros(Nin_per_f, 2, Nbin_per_x, Ntime_per_packet); % 8x2x25x20
switch data_flag
case 1 % Send white noise
sig = sqrt(sigma2/2);
......@@ -286,7 +317,15 @@ for mcnt = [0:201, 400:50:20000] %while mcnt <= 10000
data(:, 2, bin_idx, :) = CEN_imag(f_idxs,t_idxs);
data(data < 0) = 2^8 + data(data < 0);
end
case 9 % Send pulsar data
t_idxs = mod(mcnt*20 + 1:(mcnt+1)*20, Ntime);
t_idxs(t_idxs == 0) = Ntime;
f_idxs = (fid - 1)*8+1:fid*8;
freq_idxs = 5*(xid-1) + [1:5, 101:105, 201:205, 301:305, 401:405];
tmp = pulseData(f_idxs, freq_idxs, t_idxs);
data(:,1,:,:) = real(tmp);
data(:,2,:,:) = imag(tmp);
data(data < 0) = 2^8 + data(data < 0);
otherwise % Send all zeros
end
data = uint8(data);
......
fs = 155e6; % Sampling frequency - used for noise level
Ninputs = 40; % Number of inputs/antennas
Nbins = 500; % Total number of frequency bins
Nfft = 512; % F-engine FFT size
Nfengines = 5; % Number of F-engines
Nxengines = 20; % Number of X-engines (i.e. Number of GPUs)
Nin_per_f = Ninputs/Nfengines; % Number of inputs per F-engine
Nbin_per_x = Nbins/Nxengines; % Number of bins per X-engine
Ntime_per_packet = 20; % Number of time samples (spectra snapshots) per packet
Ntime = 4000;
% Increase the range of tau when dispersion measure causes m_D to exceed
% time samples.
D = 10; % DM; 10 with these parameters gives a fairly fast pulsar
freq = (0:499)*(303e3) + 1300e6; % All frequencies
fo = freq(floor(length(freq)/2)); % Center frequency
m_D = 4.1488e3*((fo^-2)-(freq.^-2))*D; % Frequency dependent timing offset
tau = -2.8e-14:((2.5e-14)+(2.8e-14))/(Ntime-1):2.5e-14; % Range of timing offsets
% tau = -3e-14:(max(m_D)+3e-14)/(Nbins-1):max(m_D); % Range of timing offsets
pulseData = zeros(Ninputs, Nbins, Ntime);
% Noisy environment
for ii = 1:size(pulseData,3)
for jj = 1:size(pulseData,2)
pulseData(:,jj,ii) = 0.1*(randn(1) + 1j*randn(1));
end
end
% for cyc = [0,1000,2000,3000] % Cycles
% for k = 1:Nbins
% tmp_tau = abs(m_D(k)-tau);
% [~, t_idx(k)] = min(tmp_tau);
% pulseData(:,k,t_idx(k)+cyc) = exp(1j*pi/4)*2/sqrt(2); % 1+1j
% pulseData(:,k,(t_idx(k)+cyc)+4) = exp(1j*pi/4); % 0.707+0.707j
% if ((t_idx(k)+cyc)-4) > 0
% pulseData(:,k,(t_idx(k)+cyc)-4) = exp(1j*pi/4); % 0.707+0.707j
% end
% end
% end
% Pulsar
pulse = 1;
for cyc = [0, -1000, 1000]
for m = 1:Ninputs
for k = 1:Nbins
tmp = abs(m_D(k)-tau);
[~,idx] = min(tmp);
phi = m*2*pi*freq(k)*tau(idx+cyc);
pulseData(m,k,idx+cyc) = pulse*exp(1j*phi) + 0.1*(randn(1) + 1j*randn(1));
end
end
end
%%%%% In switch case statement
% t_idxs = mod(mcnt*20 + 1:(mcnt+1)*20, PNtime);
% t_idxs(t_idxs == 0) = PNtime;
% f_idxs = (fid - 1)*8+1:fid*8;
% freq_idxs = 5*(xid-1) + [1:5, 101:105, 201:205, 301:305, 401:405];
% tmp = pulseData(f_idxs, freq_idxs, t_idxs);
% data(:,1,:,:) = real(tmp);
% data(:,2,:,:) = imag(tmp);
% data(data < 0) = 2^8 + data(data < 0);
%%
figure(1);
imagesc(squeeze(abs(pulseData(20,:,:))));
% imagesc(abs(exp(1j*phi)));
\ No newline at end of file
data_root = '/lustre/projects/flag';
meta_root = '/home/gbtdata';
proj_id = 'AGBT16B_400_05';
b = extract_b_output('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2017_05_03_23:17:56A.fits');
N_bin = 25;
N_pol = 4; % Order: XX*, YY*, real(XY*), and imag(XY*).
idx = 0;
for i = 1:1 % N_bin
for k = 1:1 % Order: XX*, YY*, real(XY*), and imag(XY*).
idx = idx+1;
figure(idx);
imagesc(squeeze(10*log10(abs(b(:,k,i,:)))).'); % Order: beam, pol, bin, and sti.
if k == 1
title(['X Frequency bin ' num2str(i)]);
elseif k == 2
title(['Y Frequency bin ' num2str(i)]);
elseif k == 3
title(['Real XY Frequency bin ' num2str(i)]);
elseif k == 4
title(['Imaginary XY Frequency bin ' num2str(i)]);
end
xlabel('Beam Index');
ylabel('STI Index');
end
end
\ No newline at end of file
......@@ -19,8 +19,8 @@ function [ R, dmjd, xid ] = extract_covariances( fits_filename )
Nel = 64;
Nbin = 5;
Nsamp = 40;
Nbin = 25;
Nsamp = 4000;
R = reconstruct_covariances_bdj(data, Nel, Nbin, Nsamp);
R = R(1:40, 1:40, :, :);
......
close all;
% Tone at 1367.3MHz in Bank C & D (Scan 1)
R_on_C = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2017_05_23_14:02:35J.fits');
R_on_D = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2017_05_23_14:02:35J.fits');
% No tone (Scan 2)
R_off_C = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2017_05_22_01:36:57C.fits');
R_off_D = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2017_05_22_01:36:57D.fits');
% Tone at 1367.3MHz in Bank C & D (Scan 3)
% [R_off_C,dmjd2,xid_c] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:38:36C.fits');
% [R_off_D,dmjd4,xid_d] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:38:36D.fits');
% Tone at 1367.3MHz in Bank C & D (Scan 4)
% [R_on_C,dmjd1,xid_c] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:49:39C.fits');
% [R_on_D,dmjd3,xid_d] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:49:39D.fits');
hot_c = mean(R_on_C,4);
hot_d = mean(R_on_D,4);
cold_c = mean(R_off_C,4);
cold_d = mean(R_off_D,4);
filename1 = 'Scan1_Otf_Hot_C.mat';
filename2 = 'Scan1_Otf_Hot_D.mat';
filename3 = 'Scan2_Otf_Cold_C.mat';
filename4 = 'Scan2_Otf_Cold_D.mat';
save(filename1, 'hot_c');
save(filename2, 'hot_d');
save(filename3, 'cold_c');
save(filename4, 'cold_d');
% fitsfilename = '/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:38:36D.fits';
% info = fitsinfo(fitsfilename);
% bintbl = fitsread(fitsfilename, 'binarytable', 1);
%
% data = bintbl{3};
% filename1 = 'Cold_1367MHz_D.mat';
% save(filename1, 'data');
close all; clear all;
w = zeros(40,1);
y = zeros(25,40);
% Tone at 1367.3MHz in Bank C
% [R_on_C,dmjd1,xid_c] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_14:37:11C.fits');
% [R_on_D,dmjd2,xid_d] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_14:37:11D.fits');
%
% % No tone
% [R_off_C,dmjd3,xid_c] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:01:14C.fits');
% [R_off_D,dmjd4,xid_d] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:01:14D.fits');
% % Tone at 1367.3MHz in Bank C
% % [R_on_C,dmjd1,xid_c] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:49:39C.fits');
% % [R_off_C,dmjd2,xid_c] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:38:36C.fits');
% % [R_on_D,dmjd3,xid_d] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:49:39D.fits');
% % [R_off_D,dmjd4,xid_d] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:38:36D.fits');
%
load('Scan1_Otf_Hot_C.mat');
R_on_C = hot_c;
load('Scan1_Otf_Hot_D.mat');
R_on_D = hot_d;
load('Scan2_Otf_Cold_C.mat');
R_off_C = cold_c;
load('Scan2_Otf_Cold_D.mat');
R_off_D = cold_d;
xid_c = 1;
xid_d = 2;
%
% % Bank C
% Rprime = zeros(40,40,25,1138);
% % Rprime = zeros(40,40,25);
% for i = 1:25
% for j = 1:40
% w = zeros(40,1);
% w(j) = 1;
% % hot = w(1:40)'*mean(R_on(:,:,i,:),4)*w(1:40);
% % cold = w(1:40)'*mean(R_off(:,:,i,:),4)*w(1:40);
% % hot = R_on_C(j,j,i); %mean(R_on_C(j,j,i,:),4);
% % cold = R_off_C(j,j,i); %mean(R_off_C(j,j,i,:),4);
% hot = mean(R_on_C(j,j,i,:),4);
% cold = mean(R_off_C(j,j,i,:),4);
% y(i,j) = hot/cold;
% end
% % w = ones(40,1);
% % w([3,20,30,40]) = 0;
% % w =mean(R_off(:,:,i,:),4)\w;
% % hot = w(1:40)'*mean(R_on(:,:,i,:),4)*w(1:40);
% % cold = w(1:40)'*mean(R_off(:,:,i,:),4)*w(1:40);
% % y2(i) = hot/cold;
% end
%
% for i = 1:25
% for j = 1:40
% R_on_C(3,3,:,:) = Rprime(3,3,:,:);
% R_on_C(20,20,:,:) = Rprime(20,20,:,:);
% R_on_C(30,30,:,:) = Rprime(30,30,:,:);
% R_on_C(40,40,:,:) = Rprime(40,40,:,:);
% hot_c = mean(R_on_C(j,j,i,:),4);
% cold_c = mean(R_off_C(j,j,i,:),4);
% y2(i) = hot_c/cold_c;
% % R_on_C(3,3,:) = Rprime(3,3,:);
% % R_on_C(20,20,:) = Rprime(20,20,:);
% % R_on_C(30,30,:) = Rprime(30,30,:);
% % R_on_C(40,40,:) = Rprime(40,40,:);
% % hot_c = R_on_C(j,j,i);
% % cold_c = R_off_C(j,j,i);
% y2(i) = hot_c/cold_c;
% end
% end
%
%
%
% % Bank D
% Rprime = zeros(40,40,25,1138);
% % Rprime = zeros(40,40,25);
% for i = 1:25
% for j = 1:40
% w = zeros(40,1);
% w(j) = 1;
% hot = w(1:40)'*mean(R_on_D(:,:,i,:),4)*w(1:40);
% cold = w(1:40)'*mean(R_off_D(:,:,i,:),4)*w(1:40);
% % hot_d = R_on_D(j,j,i); %mean(R_on_D(j,j,i,:),4);
% % cold_d = R_off_D(j,j,i); %mean(R_off_D(j,j,i,:),4);
% y_d(i,j) = hot/cold;
% end
% end
%
% for i = 1:25
% for j = 1:40
% R_on_D(3,3,:,:) = Rprime(3,3,:,:);
% R_on_D(20,20,:,:) = Rprime(20,20,:,:);
% R_on_D(30,30,:,:) = Rprime(30,30,:,:);
% R_on_D(40,40,:,:) = Rprime(40,40,:,:);
% hot_d = mean(R_on_D(j,j,i,:),4);
% cold_d = mean(R_off_D(j,j,i,:),4);
% y2_d(i) = hot_d/cold_d;
% % R_on_D(3,3,:) = Rprime(3,3,:);
% % R_on_D(20,20,:) = Rprime(20,20,:);
% % R_on_D(30,30,:) = Rprime(30,30,:);
% % R_on_D(40,40,:) = Rprime(40,40,:);
% % hot_d = R_on_D(j,j,i);
% % cold_d = R_off_D(j,j,i);
% % y2_d(i) = hot_d/cold_d;
% end
% end
% % % Tone at 1367.3MHz in Bank C
% % [R_on_C,dmjd1,xid_c] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:49:39C.fits');
% % [R_off_C,dmjd2,xid_c] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:38:36C.fits');
% % [R_on_D,dmjd3,xid_d] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:49:39D.fits');
% % [R_off_D,dmjd4,xid_d] = extract_covariances('/lustre/gbtdata/TGBT16A_508_01/TMP/BF/2016_10_07_15:38:36D.fits');
%
% % load('Scan1_Hot_1367MHz_C.mat');
% % R_on_C = hot_c;
% % load('Scan1_Hot_1367MHz_D.mat');
% % R_on_D = hot_d;
% % load('Scan2_Cold_NoTone_C.mat');
% % R_off_C = cold_c;
% % load('Scan2_Cold_NoTone_D.mat');
% % R_off_D = cold_d;
% % xid_c = 1;
% % xid_d = 2;
%
% % Bank C
% Rprime = zeros(40,40,25,1138);
% % Rprime = zeros(40,40,25);
% for i = 1:25
% for j = 1:40
% w = zeros(40,1);
% w(j) = 1;
% % hot = w(1:40)'*mean(R_on(:,:,i,:),4)*w(1:40);
% % cold = w(1:40)'*mean(R_off(:,:,i,:),4)*w(1:40);
% % hot = R_on_C(j,j,i); %mean(R_on_C(j,j,i,:),4);
% % cold = R_off_C(j,j,i); %mean(R_off_C(j,j,i,:),4);
% hot = mean(R_on_C(j,j,i,:),4);
% cold = mean(R_off_C(j,j,i,:),4);
% y(i,j) = hot/cold;
% end
% % w = ones(40,1);
% % w([3,20,30,40]) = 0;
% % w =mean(R_off(:,:,i,:),4)\w;
% % hot = w(1:40)'*mean(R_on(:,:,i,:),4)*w(1:40);
% % cold = w(1:40)'*mean(R_off(:,:,i,:),4)*w(1:40);
% % y2(i) = hot/cold;
% end
%
% for i = 1:25
% for j = 1:40
% R_on_C(3,3,:,:) = Rprime(3,3,:,:);
% R_on_C(20,20,:,:) = Rprime(20,20,:,:);
% R_on_C(30,30,:,:) = Rprime(30,30,:,:);
% R_on_C(40,40,:,:) = Rprime(40,40,:,:);
% hot_c = mean(R_on_C(j,j,i,:),4);
% cold_c = mean(R_off_C(j,j,i,:),4);
% y2(i) = hot_c/cold_c;
% % R_on_C(3,3,:) = Rprime(3,3,:);
% % R_on_C(20,20,:) = Rprime(20,20,:);
% % R_on_C(30,30,:) = Rprime(30,30,:);
% % R_on_C(40,40,:) = Rprime(40,40,:);
% % hot_c = R_on_C(j,j,i);
% % cold_c = R_off_C(j,j,i);
% y2(i) = hot_c/cold_c;
% end
% end
%
%
%
% % Bank D
% Rprime = zeros(40,40,25,1138);
% % Rprime = zeros(40,40,25);
% for i = 1:25
% for j = 1:40
% w = zeros(40,1);
% w(j) = 1;
% hot = w(1:40)'*mean(R_on_D(:,:,i,:),4)*w(1:40);
% cold = w(1:40)'*mean(R_off_D(:,:,i,:),4)*w(1:40);
% % hot_d = R_on_D(j,j,i); %mean(R_on_D(j,j,i,:),4);
% % cold_d = R_off_D(j,j,i); %mean(R_off_D(j,j,i,:),4);
% y_d(i,j) = hot/cold;
% end
% end
%
% for i = 1:25
% for j = 1:40
% R_on_D(3,3,:,:) = Rprime(3,3,:,:);
% R_on_D(20,20,:,:) = Rprime(20,20,:,:);
% R_on_D(30,30,:,:) = Rprime(30,30,:,:);
% R_on_D(40,40,:,:) = Rprime(40,40,:,:);
% hot_d = mean(R_on_D(j,j,i,:),4);
% cold_d = mean(R_off_D(j,j,i,:),4);
% y2_d(i) = hot_d/cold_d;
% % R_on_D(3,3,:) = Rprime(3,3,:);
% % R_on_D(20,20,:) = Rprime(20,20,:);
% % R_on_D(30,30,:) = Rprime(30,30,:);
% % R_on_D(40,40,:) = Rprime(40,40,:);
% % hot_d = R_on_D(j,j,i);
% % cold_d = R_off_D(j,j,i);
% % y2_d(i) = hot_d/cold_d;
% end
% end
%
% % figure(1);
% % plot(10*log10(abs(y2)));
% % title('Y-factor');
% % ylabel('Power');
% % xlabel('Frequency bins');
%
% figure(1);
% imagesc(10*log10(abs(y)));
% colorbar;
% title('Y-factor for XID 1');
% xlabel('Elements');
% ylabel('Frequency bins');
%
% figure(2);
% imagesc(10*log10(abs(y_d)));
% colorbar;
% title('Y-factor for XID 2');
% xlabel('Elements');
% ylabel('Frequency bins');
% % figure(1);
% % plot(10*log10(abs(y2)));
% % title('Y-factor');
% % ylabel('Power');
% % xlabel('Frequency bins');
%
% figure(1);
% imagesc(10*log10(abs(y)));
% colorbar;
% title('Y-factor for XID 1');
% xlabel('Elements');
% ylabel('Frequency bins');
%
% figure(2);
% imagesc(10*log10(abs(y_d)));
% colorbar;
% title('Y-factor for XID 2');
% xlabel('Elements');
% ylabel('Frequency bins');
figure(3);
imagesc(10*log10(abs(squeeze(mean(R_on_C(:,:,6,1),4)))));
colorbar;
title('Single frequency bin Hot-covariance');
xlabel('Elements');
ylabel('Elements');
figure(4);
imagesc(10*log10(abs(squeeze(mean(R_off_C(:,:,6,1),4)))));
colorbar;
title('Single frequency bin Cold-covariance');
xlabel('Elements');
ylabel('Elements');
bank_on = zeros(500,1);
f_idx = [1:5, 101:105, 201:205, 301:305, 401:405];
bank_on(f_idx + xid_c*5,1) = 10*log10(abs(squeeze(mean(R_on_C(4,4,:,1),4))));
bank_on(f_idx + xid_d*5,1) = 10*log10(abs(squeeze(mean(R_on_D(4,4,:,1),4))));
bank_off = zeros(500,1);
bank_off(f_idx + xid_c*5,1) = 10*log10(abs(squeeze(mean(R_off_C(4,4,:,1),4))));
bank_off(f_idx + xid_d*5,1) = 10*log10(abs(squeeze(mean(R_off_D(4,4,:,1),4))));
figure(5);
plot(bank_on, 'r');
hold on;
plot(bank_off,'b');
title('Test tone PSD');
xlabel('Frequency bins');
ylabel('Power');
legend('Hot','Cold');
hold off;
% Simple script that extracts the covariance matrices for the entire band
% and computes Y-factors
close all;
clearvars;
% System parameters
save_dir = '/lustre/projects/flag/AGBT16B_400_01/BF'; % '/lustre/projects/flag/TMP/BF'; %'/lustre/gbtdata/TGBT16A_508_01/TMP/BF';