Hi Joan,
The problem appears to be with how bodeplot is unwrapping the phase.
Here is the original code
% Delay problem
Ak = [-0.0633296088793117 + 0.00000000000000i
-0.188476918974608 + 0.00000000000000i
-0.592850411790702 + 0.00000000000000i];
Ck = [2.67546115276169e-05 + 0.00000000000000i
0.000235105989428637 + 0.00000000000000i
0.00132009482602487 + 0.00000000000000i];
D = 0;
tau = 0.001013802649678;
% Express as diagona form
F_A = eye(size(Ak,1),size(Ak,1)).*Ak;
F_B = ones(size(Ak));
F_C = Ck.';
F_D = D;
Fss = ss(F_A,F_B,F_C,F_D,'OutputDelay',tau);
% Analytic function
freq = logspace(-3,6,1000);
w = 2*pi*freq;
f = D;
for i = 1:length(Ak)
f = f + Ck(i)./(1j*w - Ak(i));
end
f = f.* exp(-1j*w.*tau);
Fana = frd(f,w);
% Comparison bode plot
opts = bodeoptions;
opts.FreqUnits = 'Hz';
opts.FreqScale = 'Log';
opts.MagUnits = 'abs';
opts.MagScale = 'Linear';
opts.grid = 'on';
figure, bodeplot(Fss,Fana,w,opts)
title("Bode plot of H")
legend('SS model', 'Analitic')
Zoom in on the end of the plot where things are different
xlim([1e4 1e6])
Now repeat the plot, but keep -180 <= phase <= 180
opts.PhaseWrapping = 'on';
figure, bodeplot(Fss,Fana,w,opts)
title("Bode plot of H")
legend('SS model', 'Analitic')
xlim([1e4 1e6])
Now both plots match.
I'm not sure why the phase unwrapping is a problem for Fana but not for Fss.
