function [an, num] = fracfilt_opt(x0, alpha, N, x_min, x_max)
  
  x = linspace(x_min, x_max, 1000);
  
  phi_approx = @(x, par) phi_hat(x, par, N);

  for i = 1:1000
    y_exact(i) = phi(x(i), x0, alpha);
  endfor
  
  % defining the initial guess
  
  for i = 1:N
    an(i) = x0 + (2 * i - 1 + alpha)/(2 * N + 1 - alpha) * (x_max - x0);
  endfor
  for i = N + 1 : 2 * N
    an(i) = x0 + (2 * (i - N) - 1 - alpha)/(2 * N + 1 - alpha) * (x_max - x0);
  endfor
  
  [y_approx, num] = leasqr(x, y_exact, an, phi_approx, 1e-4, 100);
  
endfunction


function y = phi(x, x0, alpha)

  % this function computes the value of phi
  
  y = -10 * alpha * log10(10^(2 * (x - x0)) + 1);
  
endfunction


function y = phi_hat(x, par, N)

  % this function computes the value of phi_hat
  % please notice that par(1), ..., par(N) are the zeroes, par(N + 1), ..., par(2 * N) are the poles.

  y = 0;
  for i = 1:N
    y = y + log10(10.^(2 * (x - par(i))) + 1) - log10(10.^(2 * (x - par(i + N))) + 1);
  endfor
  
  y = 10 * y;
  
endfunction
