clc; clear all;

% Define constants and parameters
MAX_TAUS = 70;
adc_clk = 133000000;    % ADC clock frequency
CIC_Rate = 64;          % CIC filter rate
dut_freq = 10000000;    % DUT frequency
adev_tap = 4;           % ADEV tap value

% Initialize Tau parameters
Tau0 = CIC_Rate / adc_clk * 10^(adev_tap + 1);

% Pre-allocate Tau structure array
Tau = repmat(struct('n', 0, 'holding', 0, 'xn', 0, 'xs', 0, 'rs', [0, 0]), MAX_TAUS, 1);

% Initialize all tau values
for i = 1:MAX_TAUS
    n = floor(10^((i+20) / 20.0)) - 9;
    Tau(i).n = n;
    Tau(i).holding = 0;
    Tau(i).xn = 0;
    Tau(i).xs = 0;
    Tau(i).rs = [0, 0];
end

% Define parameters for submitting data
% Main function to submit data
filename = 'self_corr_10M_1000s.tim';  % Replace with actual file name

% Open the file
fid = fopen(filename, 'r');
if fid == -1
    error('Failed to open file.');
end

% Skip metadata lines
while true
    line = fgetl(fid);
    if startsWith(line, 'TIC 75600 1')
        break;
    end
end

% Now, read the phase data
data_phase = [];
while true
    line = fgetl(fid);
    if ischar(line)
        % Convert the line to a number and append to data
        data_phase = [data_phase; str2double(line)];
    else
        break;  % End of file reached
    end
end
fclose(fid);

delta_phi = diff(data_phase);

%Normalize the phase data according to sampling rate
%data samples send to timelab are phase differences (/pi) between REF1 and DUT1
%calculation in Timelab is done by using sample rate & phase difference.
% pi ????????

normalized_delta_f = delta_phi*adc_clk/CIC_Rate/10^(adev_tap+1);

% Process the data
for i = 1:length(normalized_delta_f)
    y =  normalized_delta_f(i);  % Scale the frequency data
    for j = 1:MAX_TAUS
        if Tau(j).holding < Tau(j).n
            Tau(j).rs(1) = Tau(j).rs(1) + y;
            Tau(j).holding = Tau(j).holding + 1;
        elseif Tau(j).holding < 2 * Tau(j).n
            Tau(j).rs(2) = Tau(j).rs(2) + y;
            Tau(j).holding = Tau(j).holding + 1;
        else
            Tau(j).xs = Tau(j).xs + (Tau(j).rs(2) - Tau(j).rs(1))^2;
            Tau(j).rs(1) = Tau(j).rs(2);
            Tau(j).rs(2) = 0;
            Tau(j).holding = Tau(j).n;
            Tau(j).xn = Tau(j).xn + 1;
        end
    end
end

% Calculate results for each Tau
for i = 1:MAX_TAUS
    if Tau(i).xn > 0  % Check if valid
        tau(i) = Tau0 * Tau(i).n;
        var = abs(Tau(i).xs) / Tau(i).xn / Tau(i).n/Tau(i).n / 2.0;
        dev(i) = sqrt(var);
        fprintf('Tau: %.2f, Var: %.6f, Dev: %.6f\n', tau(i), var, dev(i));
    end
end

% Plot the Allan Deviation
figure(1);
loglog(tau, dev);
title('Allan Deviation of Phase Data');
xlabel('Tau (s)');
ylabel('Allan Deviation');
grid on;

N = length(data_phase);

figure(2);
plot((1:N)./20.78, data_phase,'-k');
axis tight;
grid on;
xlabel('Second  /S');
ylabel('phase difference');
ylim([min(data_phase)*1.2, max(data_phase)*1.2]);

