function [yy_sol]=Farrow_filter(M,M_high,Ns,delta,X)

% The following part calculates the Farrow coefficients from the initial
% prototype filter

N=M*Ns+M-1; % Prorotype filter order
Y=reshape(X,[M,(N+1)/2/M])';

Matr=zeros(M,M);
for br1=1:M
    for br2=1:M
        Matr(br1,br2)=(-(br2-1)/M+1/2-1/(2*M))^(br1-1);
    end
end
S=Y*(Matr)^(-1);
L=M;
M=M_high;


s0=reshape(S, [1 numel(S)]);%starting Farrow coefficients



s = sym('s',[1 numel(S)]);  % Symbolic variable s represents the Farrow coefficients that correspodents to a CMBF with a larger M
assume(s,'real')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             DEFINING MATRICES WITH COEFFICIENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ncount=M*M;
dw=0:pi/Ncount:pi/M;
Nw=length(dw);
Matr1=zeros(L,M,Nw);
for br1=1:L
    for br2=1:M
        for br3=1:Nw
        Matr1(br1,br2,br3)=(-(br2-1)/M+1/2-1/(2*M))^(br1-1)*exp(-1j*dw(br3)*(br2-1));%*z^(-m)
        end
    end
end
Matr11=zeros(L,M,Nw);%  Frequency w translated for pi/M
for br1=1:L
    for br2=1:M
        for br3=1:Nw
        Matr11(br1,br2,br3)=(-(br2-1)/M+1/2-1/(2*M))^(br1-1)*exp(-1j*(dw(br3)-pi/M)*(br2-1));%*z^(-m)
        end
    end
end

Delay2=zeros(Ns+1,L,Nw);
for br1=1:Ns+1
    for br2=1:L
        for br3=1:Nw
        Delay2(br1,br2,br3)=exp(-1j*dw(br3)*(br1-1)*M);%*z^(-nM)
        end
    end
end

Delay22=zeros(Ns+1,L,Nw);% Frequency w translated for pi/M
for br1=1:Ns+1
    for br2=1:L
        for br3=1:Nw
        Delay22(br1,br2,br3)=exp(-1j*(dw(br3)-pi/M)*(br1-1)*M);%*z^(-nM)
        end
    end
end

Ncount2=Ncount/M*10*2; % Ncount=M*M
dw2=pi/M:pi/Ncount2:pi;
Nw2=length(dw2);

Matr111=zeros(L,M,Nw2); % The matrix for the condition abs(G)<delta
for br1=1:4
    for br2=1:M
        for br3=1:Nw2
        Matr111(br1,br2,br3)=(-(br2-1)/M+1/2-1/(2*M))^(br1-1)*exp(-1j*dw2(br3)*(br2-1));%*z^(-m)
        end
    end
end

Delay222=zeros(Ns+1,L,Nw2);
for br1=1:Ns+1
    for br2=1:L
        for br3=1:Nw2
        Delay222(br1,br2,br3)=exp(-1j*dw2(br3)*(br1-1)*M);%*z^(-nM)
        end
    end
end


Matr2=zeros(L,M);
for br1=1:L
    for br2=1:M
        Matr2(br1,br2)=(-(br2-1)/M+1/2-1/(2*M))^(br1-1);
    end
end

mat_mul=zeros((Ns+1)/2,L);
for br1=1:(Ns+1)/2
    for br2=1:L
        mat_mul(br1,br2)=(-1)^(br2-1);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 END OF DEFINING MATRICES WITH COEFFICIENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

options = optimoptions('fmincon','UseParallel', true,'MaxFunEvals',1e+5,'MaxIter',1e+3,"Algorithm","sqp");
s_sol=fmincon(@(s)funkcija(s,M,L,Ns,Matr2),s0,[],[],[],[],[],[],@(s)funkcijacon(s,L,Ns,delta,Matr1,Matr11,Matr111,Delay2,Delay22,Delay222,mat_mul,Nw,Nw2),options);
s_mat_sol=reshape(s_sol,[(Ns+1)/2 L]);
g_mat_sol=s_mat_sol*Matr2;
g_sol=reshape(g_mat_sol', [1 (Ns+1)/2*M]);
gg_sol=[g_sol flip(g_sol(1:(Ns+1)/2*M))];
yy_sol=gg_sol/sum(gg_sol);% The prototype filter coefficients obtained from the Farrow coeficients

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [c]=funkcija(s,M,L,Ns,Matr2)

kk=4*M;

s_mat=reshape(s,[(Ns+1)/2 L]);
g_mat=s_mat*Matr2;
g=reshape(g_mat', [1 (Ns+1)/2*M]);
gg=[g flip(g(1:(Ns+1)/2*M))]; 
yy=gg/sum(gg); % The prototype filter coefficients obtained from the Farrow coefficients
XX1=fft(yy,kk); % Fast Fourier transform in 4M points in the frequency range [0 2*pi)

indeks2=round(kk/(M*2))+1; % Index for the stopband
indeks3=round(kk/2)+1; % Index that correspodents to pi 


fun=max(abs(XX1(indeks2:indeks3))); % 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(s,L,Ns,delta,Matr1,Matr11,Matr111,Delay2,Delay22,Delay222,mat_mul,Nw,Nw2)

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

s_mat1=reshape(s,[(Ns+1)/2 L]);



s_mat2=[s_mat1;mat_mul.*flipud(s_mat1)];% All of the coefficients of the Farrow structure


P=zeros(1,Nw*2+Nw2); % Reservation of memory for the condition functions refering to the allowed distortions of CMBF
count=1;
for br3=1:Nw
    P(count)=sum(sum((s_mat2.*Delay2(:,:,br3))*(Matr1(:,:,br3))))*conj(sum(sum((s_mat2.*Delay2(:,:,br3))*(Matr1(:,:,br3)))))+...
        sum(sum((s_mat2.*Delay22(:,:,br3))*(Matr11(:,:,br3))))*conj(sum(sum((s_mat2.*Delay22(:,:,br3))*(Matr11(:,:,br3)))))-1-delta;
    count=count+1;
    P(count)=-(sum(sum((s_mat2.*Delay2(:,:,br3))*(Matr1(:,:,br3))))*conj(sum(sum((s_mat2.*Delay2(:,:,br3))*(Matr1(:,:,br3)))))+...
        sum(sum((s_mat2.*Delay22(:,:,br3))*(Matr11(:,:,br3))))*conj(sum(sum((s_mat2.*Delay22(:,:,br3))*(Matr11(:,:,br3))))))+1-delta;
      count=count+1;
end



for br3=1:Nw2
    P(count)=abs(sum(sum((s_mat2.*Delay222(:,:,br3))*(Matr111(:,:,br3)))))-delta;
    count=count+1;
end

C=P;
Ceq=[];

end
end
