hello
try this - hope it helps
I opted for the convention of 3 dB below the peak amplitude point (which is not the 0 dB) but if you prefer -3 dB below 0 dB you can easily change that in the code at this line
threshold = 0 - 3;
instead of
threshold = max(g_dB) - 3;
full code :
%% Parameters
w = 2*pi*50; % [rad/s]
Lf = 2.5e-3; % [H]
Rf = 0.1*w*Lf; % [Ohm]
Cf = 15e-6; % [F]
Lg = 12e-3; % [H]
Rg = 0.1*w*Lg; % [Ohm]
s = tf('s');
plantg_ol = 1/((Cf*Lg)*s^2 + (Cf*Rg)*s);
plantg_cl = feedback(plantg_ol,1);
H = s*Lg + Rg;
plant_v = plantg_cl*H;
Cont_v = 0.07601*(s+2593)/s; % THIS IS GOOD!
Gv_ol = Cont_v*plant_v;
Gv_cl = feedback(Gv_ol,1);
freq = logspace(0,4,300);
[g,p] = bode(Gv_cl,2*pi*freq);
g = g(:);
g_dB = 20*log10(g);
p = p(:);
% find the BW by searching the max - 3 dB point
threshold = max(g_dB) - 3;
[ZxP,ZxN] = find_zc(freq,g_dB',threshold);
BW_3dB = ZxN(end); % -3 dB bandwith
g_3dB = interp1(freq,g_dB,BW_3dB);
p_3dB = interp1(freq,p,BW_3dB);
figure(1)
subplot(2,1,1)
semilogx(freq,g_dB,BW_3dB,g_3dB,'dr'); grid on;
title('Compensated IVL Bode Diagram');
xlabel('Frequency (Hz)')
ylabel('Gain (dB)')
l=legend([{' $G\_{c,IVL}^{CL} (s)$'};{[' BW = ' num2str(BW_3dB) ' Hz ']}]);
set(l, 'Interpreter', 'latex','FontSize',12,'Location','northeast');
% yline(0.707); xline(200);
subplot(2,1,2)
semilogx(freq,p,BW_3dB,p_3dB,'dr'); grid on;
xlabel('Frequency (Hz)')
ylabel('Phase (°)')
% yline(0.707); xline(200);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end
