% This MATLAB code plots magnitude and phase responses of the
% control-to-inductor current "transfer function" iL/ic:
% (1) exact, based on sampled-data modeling, and 
% (2) approximate, based on 1st-order or 2nd-order Pade approximations
% ECEN5807, Dragan Maksimovic, March 2007

% Buck converter example
% L = 5 uH
% m1 = m2 = 1 A/us = const.
% D = D' = 0.5
% fs = 100 kHz

fs = 1e5;
ws = 2*pi*fs;
Ts = 1/fs;
s = tf('s'); % define s variable for cont time transfer functions
z = tf('z',Ts); % define z variable for disc time transfer functions, Ts = sampling period

f=logspace(2,5,5000); % 5000 equaly spaced (log) points between 20^2=100Hz and 10^5=100 kHz
w=2*pi*f; % w is the vector of frequencies (in rad/sec)

zoh=(1-exp(-j*w*Ts))./(j*w*Ts);

% magnitude and phase responses are shown for 
% (1 - alpha)/(1-alpha e^(-sTs)) * (1-e^(-s Ts))/(s Ts)
% divising by Ts comes from sampling in front of the equivalent hold

figure(1);
clf(1,'reset');
for ramp=[0.1, 0.5, 1, 5] % ramp = ma/m2 is the parameter

    alpha = -(1-ramp)/(1+ramp);
    discrete = (1-alpha)./(1-alpha*exp(-j*w*Ts));
    H = discrete.* zoh;
    Hmagdb=20*log10(abs(H)); % compute magnitude response
    Hphu=180*unwrap(angle(H))/pi; % compute phase response

    % 1st-order approximation to e^(-sT)
    ex1 = (1 - s/(ws/pi))/(1 + s/(ws/pi));
    approx1 = ((1-alpha)/(1-alpha*ex1))*(1-ex1)/(s*Ts);
    [mag1,ph1]=bode(approx1,w);
    mag1db=20*log10(mag1);
    ph1u=unwrap(ph1);

    % 2nd-order approximation to e^(-sT)
    ex2 = (1-(pi/2)*s/(ws/2)+(s/(ws/2))^2)/(1+(pi/2)*s/(ws/2)+(s/(ws/2))^2);
    approx2 = ((1-alpha)/(1-alpha*ex2))*(1-ex2)/(s*Ts);
    [mag2,ph2]=bode(approx2,w);
    mag2db=20*log10(mag2);
    ph2u=unwrap(ph2);

    % Plot frequency responses
    subplot(2,1,1)
    %semilogx(f,Hmagdb(:,:),f,mag1db(:,:),f,mag2db(:,:)) % overlay both approx
    semilogx(f,Hmagdb(:,:),f,mag2db(:,:)) % overlay only one approx
    %semilogx(f,Hmagdb(:,:)) % H only
    grid on
    ylabel('magnitude [db]')
    axis([0 1e5 -40 20]);
    title('iL/ic magnitude and phase responses')
    hold on
    
    subplot(2,1,2)
    %semilogx(f,Hphu(:,:),f,ph1u(:,:),f,ph2u(:,:)) % overlay both approximations
    semilogx(f,Hphu(:,:),f,ph2u(:,:)) % overlay only one approx
    %semilogx(f,Hphu(:,:)) % H only
    axis([0 1e5 -180 20]);
    xlabel('frequency [Hz]')
    ylabel('phase [deg]')
    grid on
    
    hold on
end