function [a, b, c, d, A_o] = design_ellips(K, acc)
%
% function [a, b, c, d, A_o] = design_ellips(K, acc) 
% or       [num, den, A_o] = design_ellips(K, acc)   % in numerator/denominator form
% or       [num, den] = design_ellips(K, acc)        % with A_o inside numerator
%
% designs the Doppler Filter (IIR) as best as possible.
% K is the number of biquads, acc the allowed frequency error.
% Christos Komninakis, August 2001.
%
if ((nargout~=2) & (nargout~=3) & (nargout~=5)),
    disp('The number of output arguments is wrong, I cannot decide what to return!');
    return;
end;
%
inversion_flag = 0;
n = 4*K;
fD = 0.2; fDN = 2*fD;
M = 500;    % end up with M+1 frequency points, including 0 and 1.
L = floor(fDN*M);
Y = sqrt(1./sqrt(1-([0:L-1]./L).^2));
Y = [ Y sqrt(L*(pi/2-asin((L-1)/L))) ];
Y = [Y zeros(1,M-L)];
W = [0:1/M:1];
%
z = [exp(-j*W*pi); exp(-2*j*W*pi)];  % contains z_i^(-1) and z_i^(-2);
ind = [1:4:4*K]';
%
% Initialize
%
A_o = abs(randn);          % the optimum Gain coefficient.
a = randn(K,1);            % the vectors with a's, b's, c's and d's.
b = randn(K,1);
c = randn(K,1);
d = randn(K,1);
A = 100*eye(4*K);          % the matrix for the Ellipsoid Algorithm.
Upper = 1e5; Lower = 0;    % the limits for the Ellipsoid bounds.
grad_Qhat = zeros(4*K,1);
x = zeros(4*K,1);
%
% run for several outer iterations, until you get the starting point right.
%
for iteriter = 1:15,
  if (inversion_flag == 1)   % i.e if an inversion has occured in the bottom,
                             % during the previous outer iteration:
    % start ANEW, with a smaller A, and using as a starting point
    % the x you obtained after inverting the unstable zeros and poles.
    %
    inversion_flag = 0;
    A = (4*K/iteriter)*eye(4*K);   % maybe this is TOO big to restart.
    Upper = 1e3; Lower = 0;
  elseif (Upper - Lower < acc) & (Lower == 0)
    % if NO inversion occured, AND I am OK, then break.
    break;    % GET OUT, you were fine
  else
    Upper = 1e3; Lower = 0;
    A = eye(4*K);
  end;        % or do nothing, and continue Ellipsoid if you are not done yet.
  %
  % Now the inner iterations, usual Ellipsoid Method.
  %
  for iter = 1:6000,
    %
    % Compute the value of the function (Qhat) and the gradient.
    %
    H = prod(1 + [a b]*z, 1)./prod(1 + [c d]*z, 1);
    A_o = (abs(H)*Y.')/(abs(H)*abs(H.'));
    E = A_o*abs(H) - Y;
    Qhat = E*E';
    %
    grad_a = abs(H(ones(K,1),:)).*real((ones(K,1)*z(1,:))./(1+[a b]*z));
    grad_a = grad_a*E';
    grad_b = abs(H(ones(K,1),:)).*real((ones(K,1)*z(2,:))./(1+[a b]*z));
    grad_b = grad_b*E';
    grad_c = -abs(H(ones(K,1),:)).*real((ones(K,1)*z(1,:))./(1+[c d]*z));
    grad_c = grad_c*E';
    grad_d = -abs(H(ones(K,1),:)).*real((ones(K,1)*z(2,:))./(1+[c d]*z));
    grad_d = grad_d*E';
    %
    grad_Qhat(ind) = grad_a;
    grad_Qhat(1+ind) = grad_b;
    grad_Qhat(2+ind) = grad_c;
    grad_Qhat(3+ind) = grad_d;
    grad_Qhat = 2*A_o*grad_Qhat;
    %
    term = sqrt(grad_Qhat'*A*grad_Qhat);
    Upper = min(Upper, Qhat);
    Lower = max(Lower, Qhat - term);
    if (~mod(iter,200))
      disp(['Iter:', int2str(iter), ...
            '. Up. bound: ', num2str(Upper), ...
            ', Lo. bound: ', num2str(Lower), ...
            '. function Qhat: ', num2str(Qhat)]);
    end;
    %  
    if (Upper - Lower < acc)
      break;    % GET OUT of the INNER, CHECK SOLUTION STABLE.
    end;
    %
    % update
    %
    g_tilda = grad_Qhat/term;
    x(ind) = a;
    x(1+ind) = b;
    x(2+ind) = c;
    x(3+ind) = d;
    x = x - A*g_tilda/(n + 1);
    A = n^2/(n^2 - 1)*(A - 2/(n+1)*A*g_tilda*g_tilda'*A);
    a = x(ind);
    b = x(1+ind);
    c = x(2+ind);
    d = x(3+ind);
  end
  %
  % the INNER iteration is DONE here (usual Ellipsoid Algorithm).
  %
  % Now check for stability and minimum phase properties.
  %
  for i = 1:K,
    r = roots([1 a(i) b(i)]);
    viol = find(abs(r) > 1.0);
    if ( ~isempty(viol) )
      r(viol) = 1./r(viol);
      p = poly(r);
      a(i) = p(2)/p(1); b(i) = p(3)/p(1);
      inversion_flag = 1;
    end;
    %
    r = roots([1 c(i) d(i)]);
    viol = find(abs(r) > 1.0);
    if ( ~isempty(viol) )
      r(viol) = 1./r(viol);
      p = poly(r);
      c(i) = p(2)/p(1); d(i) = p(3)/p(1);
      inversion_flag = 1;
    end;
  end;
  disp(['     OUTER It:', int2str(iteriter), ...
        '. Inversion Flag:', int2str(inversion_flag)]);
  %
end  % here ends the outer iteration ("iteriter").
%
fprintf('\nDONE, Inversion Flag should be zero: %d\n', inversion_flag);
A_o = (abs(H)*Y.')/(abs(H)*abs(H.'));
E = A_o*abs(H) - Y;
Qhat = E*E';
fprintf('\nFinal value of the error, Qhat = %f\n', Qhat);
%
% Now scale for unit power:
num = [1]; den = [1];
for i=1:K,
    num = conv(num, [1 a(i) b(i)]);
    den = conv(den, [1 c(i) d(i)]);
end;
w_more = [0:0.00001:pi];
H_more = freqz(num,den,w_more);
% scale for unit power
Ao_up = A_o / sqrt(sum(abs(A_o*H_more).^2)*2*0.00001/(2*pi));
%
% Some plotting to convince everyone...
%
figure(1)
H = prod(1 + [a b]*z, 1)./prod(1 + [c d]*z, 1);
plot(W, 20*log10(A_o*abs(H)));
hold on
warning off;
plot(W, 20*log10(Y), 'r');
warning on;
hold off
grid on
title('Plot what was designed, along with the theoretical');
xlabel('Frequency, in [0, 1], where 1.0 = Nyquist Frequency');
legend('Designed', 'Theoretical', 1);
%
figure(2)
plot(W, 20*log10(Ao_up*abs(H)),'r');
hold on;
% plot it again using freqz
plot(w_more/pi, 20*log10(Ao_up*abs(H_more)),'k');
hold off;
grid on
title('Plot of what the algorithm returns, scaled for unit power');
xlabel('Frequency, in [0, 1], where 1.0 = Nyquist Frequency');
legend('Response at modeled freqs', 'Response everywhere', 1);
%
% and finally return things:
A_o = Ao_up;
if (nargout == 3),
    a = num; b = den; c = A_o;
elseif (nargout == 2),
    a = A_o*num; b = den;
end;