clear all; 
close all;

filename = "0p5.csv";     
f = 20;                   
Fs = 960;                 

data = readmatrix(filename);

% -----------------------------
% Extract channels
% -----------------------------
V = data(:,2);            % voltage column
I = data(:,3);            % current column (shunt amplifier)
N = length(V);

% -----------------------------
% Reconstruct time vector
% -----------------------------
t = (0:N-1)'/Fs;

% -----------------------------
% Remove DC offset
% -----------------------------
V_ac = V - mean(V);
I_ac = I - mean(I);

% -----------------------------
% Define sine model
% -----------------------------
w = 2*pi*f;
sine_model = @(b,t) b(1)*sin(w*t + b(2));


% Fit VOLTAGE

beta0_V = [range(V_ac)/2, 0];
beta_V = nlinfit(t, V_ac, sine_model, beta0_V);

A_out = beta_V(1);              
phi_out = beta_V(2);             % voltage phase (rad)

% -----------------------------
% Fit CURRENT
% -----------------------------
beta0_I = [range(I_ac)/2, 0];
beta_I = nlinfit(t, I_ac, sine_model, beta0_I);

A_in = beta_I(1);                % current amplitude
phi_in = beta_I(2);              % current phase (rad)

% -----------------------------
% Compute gain + phase lag
% -----------------------------
gain = A_out / A_in;             % A'/A
gain_dB = 20*log10(gain);        % Bode magnitude

phase_deg = rad2deg(phi_out - phi_in);

fprintf("Gain = %.4f (%.2f dB)\n", gain, gain_dB);
fprintf("Phase Lag = %.2f degrees\n", phase_deg);

% -----------------------------
% Plot voltage fit
% -----------------------------
t_fit = linspace(min(t), max(t), 2000);
V_fit = A_out * sin(w*t_fit + phi_out);

figure;
plot(t, V_ac, 'b'); hold on;
plot(t_fit, V_fit, 'r', 'LineWidth', 2);
xlabel("Time (s)");
ylabel("Voltage (AC)");
title(sprintf("%.2f Hz | Gain = %.3f (%.2f dB) | Phase Lag = %.2f°", ...
    f, gain, gain_dB, phase_deg));
legend("Measured Voltage", "Voltage Fit");
grid on;