1
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Just a remark concerning your last paragraph: the zero at infinity of the analog filter is not "discarded". It is mapped to the z-plane just like all the other zeros. The bilinear transform maps infinity to $z=-1$, i.e., a zero at Nyquist. $\endgroup$
    – Matt L.
    Jul 30 at 8:53
  • $\begingroup$ @MattL. I am afraid that the matlab code does indeed appear to ignore zeros at the imaginary infinite value. I updated the question with references to the lines of code that do this. $\endgroup$
    – dmedine
    Aug 2 at 0:05
  • $\begingroup$ Or, more accurately, Matlab's bilinear understands that missing zeros are infinite zeros and -1s are inserted. Which, if I am being nitpicky, is not the 'correct' implementation of a bilinear transform. It is the correct implementation of a bilinear transform when it is being used to map zeros from the s-plane to the z-plane and missing zeros are safely assumed to be infinite. $\endgroup$
    – dmedine
    Aug 2 at 0:46
1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Yes indeed bilinear.m appends -1 to the list of zeros so that its length matches the number of poles---and I see that this is the mapping from s=inf---but the result of poly is still not the expected b coefficients. $\endgroup$
    – dmedine
    Aug 2 at 0:40
  • $\begingroup$ It is for the highpass if you replace the zero at -1 with a zero at +1 $\endgroup$
    – Hilmar
    Aug 3 at 0:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.