function desGFOD(freqband,p,th,L,M,filtType,mho,pop,maxit)
%% This code is to design generalized fractional order differentiators using metaheuristic algorthm.
% p is the fractional order, it can be scalar for GFOD with fixed p or
% range vor GFOD with variable p.
% th is the range of the angle of GFOD.
% freqband is the frequency band of GFOD.
% L is the order of the subfilters.
% M is the number of sections of GFOD with variable p.
% filtType is binary to indicate whether the subfilter is IIR type of FIR type.
% mho indicates the metaheuristic algorithm used to optimize the GFOD,
% i.e., 'RCGA', 'PSO', 'WOA', or 'IWOA'.
% pop is the population size of the metaheuristic algorithm.
% maxit is the maximum number of iterations.
%% The main code
om=linspace(freqband(1)*pi,freqband(2)*pi,512);% the frequency band sample points
z=exp(1i.*om); % z variable
if length(p)==1 % check if the GFOD with fixed p
  GFODwfixedp(om,z);
else  % this for GFOD with variable p
  GFODwvarp(om,z);
end
function GFODwfixedp(om,z,th,L,filtType,mho,pop,maxit)
    F_cp=om.^p.*exp(1i*p*pi/2);% the ideal response fractional operator (positive part)
    F_cn=om.^p.*exp(-1i*p*pi/2);% the ideal response fractional operator (negative part)
    fObj=@(fn) objFn(fn,F_cp,L,z,filtType);% the objective function to be optimized
    cd (mho) % go to the folder of the optimization algorithm
    % call the function of solving this optimization problem
    if strcmp(mho,'RCGA')
        [x,err]=solveRCGA(fObj,L,pop,maxit,filtType);
    elseif strcmp(mho,'PSO')
        [x,err]=solvePSO(fObj,L,pop,maxit,filtType);
    elseif strcmp(mho,'WOA')
        [x,err]=solveWOA(fObj,L,pop,maxit,filtType);
    else
        [x,err]=solveIWOA(fObj,L,pop,maxit,filtType);
    end
    cd ('..') % return to the main folder
    % plot the figures for the new design
    c1=sin(pi/2*(p+th))/sin(p*pi); % the coefficient c1
    c2=sin(pi/2*(p-th))/sin(p*pi); % the coefficient c2
    H_d=c1'*F_cp+c2'*F_cn; % the ideal reponse of fractional operator
    H_1=c1'*Fzp+c2'*Fzn; % the transfer function of GFOD with fixed p
    subplot(3,1,2)
    plot(om,angle(Fzp)); % plot the angle of the positive part of GFOD
    hold on
    plot(om,angle(H_1(20,:))); % plot the angle of the GFOD
    % compute the error vs. theta, i.e., E1
    for k=1:length(th)
       Errth(k)=sum(abs((H_d(k,:))-(H_1(k,:).*exp(0i.*om*pi/2))).^2/512).^0.5;
    end
    subplot(3,1,3)
    plot(th,Errth) % plot the error
    
function GFODwvarp(om,z,th,L,M,filtType,mho,pop,maxit)
    % construct the ideal response ar each p and omega
    F_cp(1,:)=om.^p(1).*exp(1i*p(1)*pi/2); % the positive part of the ideal response
    F_cn(1,:)=om.^p(1).*exp(-1i*p(1)*pi/2); % the negative part of the ideal response
    for l=2:length(p)
      F_cp(l,:)=om.^p(l).*exp(1i*p(l)*pi/2); 
      F_cn(l,:)=om.^p(l).*exp(-1i*p(l)*pi/2);
    end
    fObj=@(an) objFn2(an,F_cp,M,L,z,p,filtType); % the objective function to be optimized
    cd (mho) % go to the folder of the optimization algorithm
    % call the function of solving this optimization problem
    if strcmp(mho,'RCGA')
        [x,err]=solveRCGA(fObj,(L+1)*(M+1),pop,maxit,filtType);
    elseif strcmp(mho,'PSO')
        [x,err]=solvePSO(fObj,(L+1)*(M+1),pop,maxit,filtType);
    elseif strcmp(mho,'WOA')
        [x,err]=solveWOA(fObj,(L+1)*(M+1),pop,maxit,filtType);
    else
        [x,err]=solveIWOA(fObj,(L+1)*(M+1),pop,maxit,filtType);
    end
    cd ('..') % return to the main folder
    % plot the figures
    Fzpp=trnsFn2(x,M,L,z,p,-1,filtType);% the positive part of the new GFOD
    Fzpn=trnsFn2(x,M,L,z,p,1,filtType); % the negative part of the new GFOD
    % compute the coefficients c1 and c2 for each p and theta
    c1=zeros(length(th),length(p)); 
    c2=zeros(length(th),length(p));
    for m=1:length(th)
      c1(m,:)=sin(pi/2*(p+th(m)))./sin(p*pi);
      c2(m,:)=sin(pi/2*(p-th(m)))./sin(p*pi);
    end
    % construct the ideal response and the proposed GFOD
    H_d=zeros(length(th),length(p),length(om));
    H_2=zeros(length(th),length(p),length(om));
    for l=1:length(p)
     for m=1:length(th)
        H_d(m,l,:)=c1(m,l).*F_cp(l,:)+c2(m,l).*F_cn(l,:);
        H_2(m,l,:)=c1(m,l).*Fzpn(l,:)+c2(m,l).*Fzpp(l,:);
        err2(m,l)=sum(abs(H_d(m,l,:)-H_2(m,l,:)).^2/length(om)).^0.5;
     end
        err3(l,:)=abs(F_cn(l,:)-Fzpn(l,:)).^2/length(om);
    end   
    figure
    ideal1=H_d(2,1,:);
    ideal=ideal1(:);
    plot(om,abs(ideal));
    hold on
    gfod1=H_2(2,1,:);
    gfod=gfod1(:);
    plot(om,abs(gfod));
    fig1=figure('Units','inches','Position',[0 0 4.08 1.92],'PaperPositionMode','auto');
    axes1 = axes('Parent',fig1);
    hold(axes1,'on');
    load('myCustomColormap.mat')
    [P,Th]=meshgrid(p,th);
    surf(P,Th,err2,'Parent',axes1,'LineWidth',1);
    colormap(myCustomColormap);
    xlabel('$p$','interpreter','latex');
    ylabel('$\theta$','interpreter','latex');
    zlabel('Error');
    xlim([0 1]);
    ylim([-2 2]);
    view(axes1,[-145 35]);
    box(axes1,'on');
    grid(axes1,'on');
    % Set the remaining axes properties
    set(axes1,'FontName','Times','FontSize',8,'GridAlpha',0.9,'GridLineStyle',...
    ':');
function err=objFn(fn,F_c,L,z,FIR)
% objective function for GFOD with fixed p
  Fz=trnsFn(fn,L,z,-1,FIR);
  err=sum((abs((F_c)-(Fz)).^2)/512).^0.5;

function Fz=trnsFn(fn,L,z,sign,FIR)
% compute the transfer function for the GFOD with fixed p
if(nargin<5)
    FIR=0;
end
Fz=0;
Bz=0;
Az=0;
if (FIR==1)
for n=-L:L
    Fz=Fz+fn(n+L+1).*z.^(sign*n);
end
else
for n=0:L
    Bz=Bz+fn(n+1).*z.^(sign*n);
end
for n=0:L
    Az=Az+fn(n+L+2).*z.^(sign*n);
end
Fz=Bz./Az;
end
function err=objFn2(an,F_c,M,L,z,p,FIR)
% objective function for GFOD with variable p
Fzp=trnsFn2(an,M,L,z,p,1,FIR);
err=sum(sum((abs(F_c-Fzp)).^2/512)/9).^0.5;

function Fzp=trnsFn2(an,M,L,z,p,sign,FIR)
% transfer function of GFOD with variable p
if(nargin<5)
    FIR=0;
end
Fzp=0;
if (FIR==1)
    for k=0:M
      Az=0;
      for n=-L:L
        Az=Az+an(n+L+1+2*L*k+k).*z.^(sign*n);
      end
      Fzp=Fzp+(p.^k)'*Az;
    end
else
  for k=0:M
      Bz=0;
      Az=0;
      for n=0:L
         Bz=Bz+an(n+1+2*k*L+2*k).*z.^(sign*n);
      end
      for n=0:L
         Az=Az+an(n+L+2+2*k*L+2*k).*z.^(sign*n);
      end
      Fzp=Fzp+(p.^k)'*(Bz./Az);
  end
end