%% example 1
    clear
    close all
    clc

    Ts = 1;           % sampling period
    alpha = -0.5;     % fractional order

    om = logspace(-4,pi,451);
    t_max = 31000;

%% ideal characteristic of FO integrator
    H(1,1,:) = (1 - exp(-1i*om)).^alpha;
    imp_idealI(:,1) = cumsum([1 -parPJ(-alpha,t_max+1)]);


%% integer-order approximation of a fractional-order integrator via Takenaka-Malmquist functions

    w_L = 1.3237e-05;
    w_H = 6.3427;       % w_L and w_H - lower and higher frequency bounds of approximation
    M = 11;             % number of filters
   
    [A,B,C,D] = FODifTM(alpha+1,M,w_L,w_H,Ts);
    sys = ss(A,B,C,D,Ts);
    int = ss(1,1,1,1,Ts);
    sysI = series(sys,int);
    Htm = freqresp(sysI,om);

    
%% Discrete Oustaloup approximation of a s^r operator   
% required FOMCON toolbox for MATLAB (Aleksei Tepljakov)
% https://www.mathworks.com/matlabcentral/fileexchange/66323-fomcon-toolbox-for-matlab
% modification: Check upper frequency bound disabled  
        
    w_L = 8.2683e-06;
    w_H = 5.3040;         % w_L and w_H - lower and higher frequency bounds of approximation

    N = floor(M/2);     % N relates to the number of  poles: M = 2N+1,       
    oust = ss(doustafod(alpha+1,N,w_L,w_H,Ts));
    oustI = series(oust,int);
    Houst = freqresp(oustI,om);
    
    freq_err_TM = sqrt(sum(abs(H(:) - Htm(:)).^2))

    freq_err_O  = sqrt(sum(abs(H(:) - Houst(:)).^2))
    
    
    
    
    %% plots
    clear om H 
    
    om = logspace(-7,pi,751);
    H = (1 - exp(-1i*om)).^alpha;
  
    [magtm,phtm] = bode(om,sysI);
    [mago,pho] = bode(om,oustI);
    
    subplot(2,1,1)
    semilogx (om,20*log10(magtm(:)),'b');
    hold on; grid;  x=axis;x(2)=pi;axis(x);
    semilogx (om,20*log10(mago(:)),'r');
    semilogx (om,20*log10(abs(H)),'g');
    legend("T-M","Oustaloup","Actual",'Location','southeast','Interpreter','latex','FontSize',9);
    x=axis;x(2)=pi;axis(x);
    title("M = " + M +", $\alpha$ = " + alpha,'interpreter','latex');
        
    subplot(2,1,2)
    semilogx (om,phtm(:),'b');
    hold on; grid;  x=axis;x(2)=pi;axis(x);
    semilogx (om,pho(:),'r');
    semilogx (om,180/pi*unwrap(angle(H)),'g');
    
    
    
    
    imp_tmI   = impulse(sysI,t_max);
    imp_oustI = impulse(oustI,t_max);
    
    time_err_TM = sqrt(sum((imp_idealI - imp_tmI).^2))
    time_err_O  = sqrt(sum((imp_idealI - imp_oustI).^2))
    
    
    