how to convert the zeros of an odd order Chebyshev Type II from low-pass prototype?

Question

My DSP book (Rabiner, Lawrence R., and Bernard Gold. "Theory and application of digital signal processing." Englewood Cliffs: Prentice-Hall (1975).) has the following transformations for mapping a lowpass (analog) prototype filter to the desired filter:

Lowpass to lowpass: $$ s_p\rightarrow\frac{s}{\Omega_u} $$ Lowpass to highpass: $$ s_p\rightarrow\frac{\Omega_u}{s} $$ Lowpass to bandpass: $$ s_p\rightarrow\frac{s^2 + \Omega_l\Omega_u}{s(\Omega_u-\Omega_l)} $$ Lowpass to notch: $$ s_p\rightarrow\frac{s(\Omega_u-\Omega_l)}{s^2 + \Omega_l\Omega_u} $$ where $s_p$ is the poles and zeros of the prototype, $s$ is the poles and zeros of the desired filter and $\Omega$ are the desired cutoff frequencies.

I want to compute filter coefficients for a digital IIR filter using the Chebyshev Type II design technique in a C++ library and am verifying my results against Matlab. I have methods for calculating the prototype, computing the transformations, doing the bilinear transform and multiplying out the polynomials. It works perfectly so long as the number of poles and zeros is even.

I am 'checking by hand' with the following Matlab code:

Wo = .5; % cutoff as percentage of nyquist frequency
fs = 2; % matlab convention for working with digital frequencies and analog filter prototypes
uo = 2*fs*tan(pi*Wo/fs); 
order = 4;

[z,p,k] = cheb2ap(order,40); % compute prototype

zt = uo*z; % mapping low pass prototype to low pass at desired cutoff
pt = uo*p;

[ztt,ptt,ktt] = bilinear(zt,pt,k,fs); % analog to digital

btt = poly(ztt)*ktt; % multiply the roots of the z-plane polynomial to get the filter coefficients
att = poly(ptt);

[b,a] = cheby2(order,40,Wo); % check against the built in Matlab function

This yields identical values for the feed forward filter coefficients. The zeros are also (obviously) the same but in a different order.

>> b/btt
ans =

  1.000000000000001 - 0.000000000000000i

And such is the case when the order (I know it isn't the actual order of the filter, it is the order - 1) is any even value. However, if I set order to any odd value, the result is very nearly 4.0, which is the value I get for $\Omega_u$.

Ok, so I don't understand that at all, but at least I can work with it.

However, if I instead map to another filter type, e.g. highpass, I also get the same poles and zeros for even numbers for order, but for odd numbers, the zeros are off by a different value whose meaning I can't fathom (poles are always the same):

Wo = .5;
fs = 2;
uo = 2*fs*tan(pi*Wo/fs);
order = 3;

[z,p,k] = cheb2ap(order,40);

zt = uo./z; % highpass mapping
pt = uo./p;

[ztt,ptt,ktt] = bilinear(zt,pt,k,fs);

btt = poly(ztt)*ktt;
att = poly(ptt);

[b,a] = cheby2(order,40,Wo, 'high')

b/btt

ans =

     2.440041381447946e-13

Similarly with the bandpass and notch cases (which require solving a quadratic to map the poles and zeros---so I leave it out for simplicity's sake). What's worse is that in this case, the quotient of b/btt changes depending on the value of order:

...
order=5;
...
b/btt
ans =

     -1.672191186690123e-13 - 1.165049447799913e-29i

...
order=7;
...
b/btt
ans =

    2.037229248854087e-13 + 1.050686044147603e-29i

We all know that in Chebyshev filters, if the number of zeros is odd, there will be a zero at infinite where the denominator is 0. In the Matlab's bilinear transform, this imaginary, infinite, zero gets discarded by the bilinear transform function (lines 101-102 of bilinear.m):

zs = zs(isfinite(zs));  % Strip infinities from zeros
                            % Do bilinear transformation

but this doesn't matter because cheb2ap skips the infinite zeros anyway (lines 29-35 of cheb2ap):

if isodd(n)
    m = n - 1;
    z = cos([1:2:n-2 n+2:2:2*n-1]*pi/(2*n))';
else
    m = n;
    z = cos((1:2:2*n-1)*pi/(2*n))';
end

and at the end of bilinear.m (line 118), a zero of value -1 is appended to the list of zeros so that the length matches the list of poles:

zd = [zd1;-ones(length(pd)-length(zd1),1)];

I suspect that here is where I need to make some kind of adjustment, but I don't know what it should be. Any insight is much appreciated.

Answer 1

This is a bit easier to see if you look at the zeros directly and not at the numerator polynomial, i.e. [z1,p1,k1] = cheby2(order,40,Wo, 'high');

Odd order lowpass filters have a zero at $s=\infty$ or $z=-1$. The prototype function does NOT include this zero. For example buttap does not return any zeros at all.

Apparently bilinear() notices the missing zeros and adds them back assuming that they were are $z=-1$.

For a high pass filter you have to manually flip the sign of that zero.