clear all;

fid = fopen('coeffReal.dat');

[realCoeff, count] = fscanf(fid, '%f');

fclose(fid);

fid = fopen('coeffImag.dat');

[imagCoeff, count] = fscanf(fid, '%f');

fclose(fid);

branchNo = 13;
sampleNo = 5000;

data = complex(realCoeff, imagCoeff);
clear readCoeff imagCoeff;

data = reshape(data, [branchNo, sampleNo]).';


fprintf(1, '\n');

dc = mean(data);
for ii=1:branchNo
  %  fprintf(1, 'DC(%d) is %f+j*%f\n', ii, real(dc(ii)), imag(dc(ii)));
end

fprintf(1, 'The difference between DC and 0 is %f.\n\n', norm(dc));

% fprintf(1, 'Correlation Matrix:');
corr = data'*data/sampleNo;

fprintf(1, 'The difference between measured spatial correlation matrix and I is %f\n', norm(corr-eye(branchNo)));


clf;
windowSize = 17;
aux = data(:, 5);
aux = [aux(2:end); aux(1)]; % rotate aux by one
selectedBranchIndex = 10;
temporalCorrelation = xcorr(data(:, selectedBranchIndex), data(:, selectedBranchIndex), windowSize, 'coeff');
mag = abs(temporalCorrelation);
plot(mag);

fd = 100*18.46/1000/4;
theoretical = besselj(0, 2*pi*fd*(-windowSize:windowSize));
hold on;
plot(abs(theoretical), 'r');
legend('measured correlation', 'theoretical correlation');

fprintf(1, '\nSelect branch # %d to calculate the temporal correlation.\n', selectedBranchIndex);
fprintf(1, 'The difference between measured temporal correlation and Bessel function is %f within window size of %d.\n', norm(temporalCorrelation-theoretical.'), windowSize);