function [A0d,A0n,A1d,A1n] = emqf_des(N,alpha,Fp,Fs)
% This program designs EMQF filter for the following input parameters:
% N - filter order
% alpha - common constants in the second order sections
% Fp - Passband cut-off frequency
% Fs - Stopband cut-off frequency
% Program returns the coefficients of the all-pass branches A0(z) and A1(z)
%
% Designing the start-up halfband filter
ksi = tan(pi*Fs/2)/tan(pi*Fp/2); % selectivity factor
Fp_HB = 2*atan(1/sqrt(ksi))/pi; % passband edge frequency of the halfband filter 
[b,a,z,p,k] = halfbandiir(N,Fp_HB); % halfband filter design

pp = sort(abs(p).^2);

% Computing the coefficients of halfband filter all-pass branches
alphaHB = 0;
alpha1HB = 0;
beta_HB = pp(2:2:length(pp));

% EMQF filter

% Frequency transformations
% alpha = -(1-1/8-1/64);
if alpha == 0, alpha1=0; else
alpha1 = (1-sqrt(1-alpha^2))/alpha;
end
beta = (beta_HB+alpha1^2)./(beta_HB*alpha1^2+1);

% Computing the coefficients of EMQF filter all-pass branches
p1 = [1 alpha1];
for ind = 2:2:(N-1)/2
  p1 = conv(p1,[1 alpha*(1+beta(ind)) beta(ind)]);
end
p2 = 1;
for ind = 1:2:(N-1)/2
  p2 = conv(p2,[1 alpha*(1+beta(ind)) beta(ind)]);
end

% EMQF filter all-pass branches
A1d = p1;
A1n = fliplr(A1d);
A0d = p2;
A0n = fliplr(A0d);