
function [X]=Function_evenM_oddN(M,N,delta1,delta2,delta3,ratio)

Fpass = 0.6*1/M;  % Passband frequency of the prototype filter
Fstop = 1/M;  % Stopband frequency of the prototype filter
Dpass = 1;   % Passband ripple for the design od of the optimisation starting filter
Dstop= 1*ratio;   % Stopband attenuation for the design od of the optimisation starting filter
dens  = 20;    % Density factor

% Calculate the coefficients using the FIRPM function.
[N1, Fo, Ao, W] = firpmord([Fpass, Fstop], [1 0], [Dpass, Dstop]);
b  = firpm(N, Fo, Ao, W, {dens}); % Note that filter order is N, not N1
x = sym('x',[1 (N+1)/2]); % Symbolic variable x represents the first (N+1)/2 coefficients of the prototype filter while it has symetric impulse response, 
                          % so only (N+1)/2 coefficients are unknown. 



assume(x,'real')  % Variable x is real


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% DEFINING MATRICES WITH COEFFICIENTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha0=exp(1i*pi/4);
beta0=exp(-1i*(N)*pi/(4*M));
NN=N+1; % The number of the eqidistance disctere frequencies in the range [0 pi]
niz1=zeros((N+1)/2,M,M,NN+1); % Reservation of memory
niz2=zeros((N+1)/2,M,M,NN+1); % Reservation of memory
niz3=zeros((N+1)/2,M,M,NN+1); % Reservation of memory
niz4=zeros((N+1)/2,M,M,NN+1); % Reservation of memory
dw=pi/NN;

% The following part calculates the coefficients that form the distortion transfer function and the
% aliasing transfer function.
% br1 is a time instance, br2 referrs to a chanel (1. to Mth channel), br3
% inexes is the condition function in the
% constrained optimisation and br4 refers to the disctere frequencies.


for br1=1:(N+1)/2 
    
    for br2=1:M 
        
        for br3=1:M
           for br4=1:NN+1
               if br3==1

                             niz1(br1,br2,br3,br4)=(exp(-1*1i*(br1-1)*((br4-1)*dw-(2*br2-1)*pi/(2*M))) ...
                                 +exp(-1*1i*(N-(br1-1))*((br4-1)*dw-(2*br2-1)*pi/(2*M))))*alpha0^((-1)^(br2-1))*beta0^(2*(br2-1)+1);
                             niz2(br1,br2,br3,br4)=(exp(-1*1i*(br1-1)*((br4-1)*dw+(2*br2-1)*pi/(2*M))) ...
                                 +exp(-1*1i*(N-(br1-1))*((br4-1)*dw+(2*br2-1)*pi/(2*M))))*conj(alpha0^((-1)^(br2-1)))*conj(beta0^(2*(br2-1)+1));
                             niz3(br1,br2,br3,br4)=(exp(-1*1i*(br1-1)*((br4-1)*dw-(2*br2-1)*pi/(2*M))) ...
                                 +exp(-1*1i*(N-(br1-1))*((br4-1)*dw-(2*br2-1)*pi/(2*M))))*conj(alpha0^((-1)^(br2-1)))*beta0^(2*(br2-1)+1);
                             niz4(br1,br2,br3,br4)=(exp(-1*1i*(br1-1)*((br4-1)*dw+(2*br2-1)*pi/(2*M))) ...
                                 +exp(-1*1i*(N-(br1-1))*((br4-1)*dw+(2*br2-1)*pi/(2*M))))*alpha0^((-1)^(br2-1))*conj(beta0^(2*(br2-1)+1));

               else
                   l=br3-1;

                             niz1(br1,br2,br3,br4)=(exp(-1*1i*(br1-1)*((br4-1)*dw+(-2*pi*l/M)-(2*br2-1)*pi/(2*M))) ...
                                 +exp(-1*1i*(N-(br1-1))*((br4-1)*dw+(-2*pi*l/M)-(2*br2-1)*pi/(2*M))))*alpha0^((-1)^(br2-1))*beta0^(2*(br2-1)+1);
                             niz2(br1,br2,br3,br4)=(exp(-1*1i*(br1-1)*((br4-1)*dw+(-2*pi*l/M)+(2*br2-1)*pi/(2*M))) ...
                                 +exp(-1*1i*(N-(br1-1))*((br4-1)*dw+(-2*pi*l/M)+(2*br2-1)*pi/(2*M))))*conj(alpha0^((-1)^(br2-1)))*conj(beta0^(2*(br2-1)+1));
                             niz3(br1,br2,br3,br4)=(exp(-1*1i*(br1-1)*((br4-1)*dw-(2*br2-1)*pi/(2*M))) ...
                                 +exp(-1*1i*(N-(br1-1))*((br4-1)*dw-(2*br2-1)*pi/(2*M))))*conj(alpha0^((-1)^(br2-1)))*beta0^(2*(br2-1)+1);
                             niz4(br1,br2,br3,br4)=(exp(-1*1i*(br1-1)*((br4-1)*dw+(2*br2-1)*pi/(2*M))) ...
                                 +exp(-1*1i*(N-(br1-1))*((br4-1)*dw+(2*br2-1)*pi/(2*M))))*alpha0^((-1)^(br2-1))*conj(beta0^(2*(br2-1)+1));

               end

           end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X0=b(1:(N+1)/2);

options = optimoptions('fmincon','UseParallel', true,'MaxFunEvals',1e+5,"Algorithm","sqp");

X=fmincon(@(x)funkcija(x,N,M),X0,[],[],[],[],[],[],@(x)funkcijacon(x,N,M,delta1,delta2,delta3,niz1,niz2,niz3,niz4),options); 
X=[X flip(X(1:(N+1)/2))];
X=X/sum(X); % Coefficients of the prototype filter
% First (N+1)/2 coefficients of the prototype filter which are the result
% of the constraint optimisation. funkcija(x,N,M) is the goal function, and
% funkcijacon is the constraint function

function [c]=funkcija(x,N,M)
xx=[x flip(x(1:(N+1)/2))]; % In xx are all the prototype filter coefficients
yy=xx/sum(xx); % Normalisation of xx
XX1=fft(yy,(N+1)*4); % Frequency resolution is N*4, XX1 is the Fast Fourier transform of the normalised filter coefficients
index2=round(4*(N+1)/(2*M))+1; % Index for the stopband
index3=round(4*(N+1)/2)+1; % Index that correspodents to pi 

fun=max(abs(XX1(index2:index3))); % Maximum of the amplitude spectrum of the prototype filter in she stopband region

c=fun; % Output value of the goal function

end


function [C,Ceq]=funkcijacon(x,N,M,delta1,delta2,delta3,niz1,niz2,niz3,niz4)


NN=N+1; % The number of the eqidistance disctere frequencies in the range [0 pi]

xx=[x flip(x(1:(N+1)/2))]; % All coeficients of the prototype filter
yy=xx/sum(xx); % Normalisation of the coefficients
XX1=fft(yy,N+1); % Fast Fourier transform in N points in the frequency range [0 2*pi)


index2=round((N+1)/(2*M))+1; % Index that correspodents to the stopband of the prototype filter
index3=round((N+1)/2)+1; % Index that correspodents to pi 

pomocna=zeros(1,(NN+1)*(M+1)+1+index3-index2); % Reservation of memory

brojac=1;
for br=1:M 
    if br==1
        for brw=1:(NN+1)
         % S1 and S2 are the condition functions refering to the distortion transfer function
         S1=(abs((1/M)*sum((x(1:((N+1)/2))*niz3(:,:,br,brw)+x(1:(N+1)/2)*niz4(:,:,br,brw)).*(x(1:(N+1)/2)*niz1(:,:,br,brw)+x(1:(N+1)/2)*niz2(:,:,br,brw))))-1-delta1);
         S2=(-1*abs((1/M)*sum((x(1:(N+1)/2)*niz3(:,:,br,brw)+x(1:(N+1)/2)*niz4(:,:,br,brw)).*(x(1:(N+1)/2)*niz1(:,:,br,brw)+x(1:(N+1)/2)*niz2(:,:,br,brw))))+1-delta1);
         pomocna(brojac)=S1;
         pomocna(brojac+1)=S2; 
         brojac=brojac+2; 
        end
    else
        for brw=1:(NN+1)
        % S is the condition function refering to the aliasing transfer function
        S=(abs((1/M)*sum((x(1:(N+1)/2)*niz3(:,:,br,brw)+x(1:(N+1)/2)*niz4(:,:,br,brw)).*(x(1:(N+1)/2)*niz1(:,:,br,brw)+x(1:(N+1)/2)*niz2(:,:,br,brw))))-delta2);
        pomocna(brojac)=S;
        brojac=brojac+1;

        end
    end
end


for brr=brojac:brojac+index3-(index2) % brr coresponds to the stopband frequencies of the prototype filter
pomocna(brr)=abs(XX1(index2+brr-brojac))-delta3; % Condition function considering stopband attenuation of the prototype filter
end

C=pomocna; % Output of the optimisation constraint function 
Ceq=[];

end
end
