clear all; 
close all;

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

data = readmatrix(filename);

V = data(:,2);            
I = data(:,3);            
N = length(V);

t = (0:N-1)'/Fs;
% remove offset for values
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);             

% 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);

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;
