Phase changing Allpass IIR Filter

The following Matlab function generates allpass coefficients for an IIR filter.  In this design, the magnitude response is unchanged but the phase response is very different.  This code only supports 1st order or 2nd order variants.

1st order allpass filters shift the phase by 180 degrees, while the 2nd order shifts the phase by 360 degrees.  The "center frequency" of Fc defines where the phase response should be halfway to the max shift.

For example,

If order=2 and Fc = 2000Hz, there would be a 180deg shift at 2kHz

If order=1 and Fc = 5000Hz, there would be a 90deg shift at 5kHz

function [b,a] = allpass(order, Fc, Fs, Q)
% Returns allpass filter coefficients.
% Currently only 1st and 2nd orders are supported.
%
% Usage: [b,a] = ALLPASS(N, FC, FS, Q);
%
%        N is the order of the allpass
%        FC is the frequency a the 90deg phase shift point
%        FS is the sampling rate
%        Q is quality factor describing the slope of phase shift
%
% Author: sparafucile17 01/07/2004

if(order > 2)
    error('Only 1st and 2nd orders are supported');
end

%Bilinear transform
g = tan(pi*(Fc/Fs));

if(order == 2)
    d  = 1/Q;
    K  = 1/(1 + d*g + g^2);
    b0 = (1 - g*d + g^2) * K;
    b1 = 2 * (g^2 - 1) * K;
    b2 = 1;
    a1 = b1;
    a2 = b0;
else
    b0 = (g-1)/(g+1);
    b1 = 1;
    b2 = 0;
    a1 = b0;
    a2 = 0;
end

b = [b0 b1 b2];
a = [1  a1 a2];