IIR allpass filter design using an equation error method

iir_ap.m

function a = iir_ap(N,f,phi,W)
% function a = iir_ap(N,f,phi,W)
% IIR allpass filter design using an equation error method
%
% a		denominator coefficients; numerator = flipud(a)
% N		filter order
% f		frequency grid; 0 <= f <= 1 ('1' equals Nyquist)
% phi	desired phase on the grid (in radians)
% W		positive weight function on the frequency grid
%
% M. Lang, mattsdspblog@gmail.com

% EXAMPLE 1: fractional delay of 7.5 samples
% f = linspace(0,.95,200);
% w = pi*f;
% phi = -7.5*w;
% W = ones(1,200);
% N = 8;
% a = iir_ap(N,f,phi,W);

% EXAMPLE 2: equalization of phase of Cauer lowpass filter
% [bc,ac] = ellip(5,.5,50,.5);
% [Hc,w] = freqz(bc,ac,128);
% W = ones(128,1); W(1:65) = 10;
% Pc = angle(Hc);
% del = 10;
% phi = - Pc - del * w;
% N = 12;
% a = iir_ap(N,w/pi,phi,W);

% EXAMPLE 3: Hilbert transformer
% f = linspace(0,1,200);
% w = pi*f;
% W = ones(1,200);
% del = 9;
% phi = - pi/2 - del * w;
% N = 10;
% a = iir_ap(N,f,phi,W);

w = pi*f(:);
phi = phi(:);
W = sqrt( W(:) );
A = exp( - 1i * w * ( N-1 : -1 : 0 ) ) - exp( 1i * ( phi(:,ones(1,N)) - w * (1:N) ) );
b = exp( 1i * phi ) - exp( - 1i * N * w );
A = W(:, ones(1,N) ) .* A;
b = W .* b;
a = [ real(A); imag(A) ] \ [ real(b); imag(b) ];
a = [ 1; a ];