I am currently trying to implement a C++ program to calculate digital filter coefficients using a Chebyshev Type II filter design technique. The text I am using to understand the theory of the technique is chapter 4.9, section 3 from Rabiner and Gold (Rabiner, Lawrence R., and Bernard Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs: Prentice-Hall, 1975, pp230-235). I am comparing their approach with Matlab's cheby2 function which comes from Parks, Thomas W., and C. Sidney Burrus. Digital Filter Design. New York: John Wiley & Sons, 1987, chap. 7. (which I don't have). When calculating the analog prototype filter (cheb2ap in Matlab) I find that Rabiner and Gold define the same zeros but different poles. I am looking for an explanation of what this difference is the result of and what the theoretical implications are down the line (i.e. transforming the prototype to a new analog filter, digitizing, and converting the roots into polynomials to get the final result) as I cannot reason it out myself.

I include the equations from Rabiner and Gold as well as the Matlab code that computes the poles as well as a transcription into mathematical formulae (which I hope I got right).

Rabiner and Gold state on p233 that the poles in Chebyshev Type II filters are given by: $$ s_k = \sigma_k + j\Omega_k, k= 1,2,\ldots,n $$

where $$ \sigma_k = \frac{\Omega_r\alpha_k}{\alpha^2_k + \beta^2_k}\\ \Omega_k = \frac{-\Omega_r\beta_k}{\alpha^2_k + \beta^2_k} $$ and $$ \alpha_k = -\sinh\varphi\sin[\frac{(2k-1)\pi}{2n}]\\ \beta_k = \cosh\varphi\cos[\frac{(2k-1)\pi}{2n}] $$ and $$ \sinh\varphi = \frac{\gamma - \gamma^{-1}}{2}\\ \cosh\varphi = \frac{\gamma + \gamma^{-1}}{2} $$ and $$ \gamma = \left(A + \sqrt{A^2-1}\right)^{1/n} $$ My parameters are $A=40$ (this is the stopband loss parameter in dB), $n=6$ and since this is my prototype $\Omega_r = 1$.

The poles yielded by these equations is $$ s_1 = -0.131536256147682 - j*0.787637781759915\\ s_2 = -0.496441993970881 - j*0.796529798508779\\ s_3 = -1.096353099998444 - j*0.471342301058600\\ s_4 = -1.096353099998444 + j*0.471342301058600\\ s_5 = -0.496441993970882 + j*0.796529798508779\\ s_6 = -0.131536256147682 + j*0.787637781759915 $$

The Matlab code for finding the poles is the following (from cheb2ap):

delta = 1/sqrt(10^(.1*rs)-1);
mu = asinh(1/delta)/n;
p = exp(1i*(pi*(1:2:2*n-1)/(2*n) + pi/2)).';
realp = real(p); realp = (realp + flipud(realp))./2;
imagp = imag(p); imagp = (imagp - flipud(imagp))./2;
p = complex(sinh(mu)*realp, cosh(mu)*imagp);

where rs is $A$ and n is $n$. The cutoff frequency is implicitly set to 1 rad/second. Or in formulae: $$ \mu = \frac{\sinh^{-1}(1/\delta)}{n} $$ where $$ \delta = \frac{1}{\sqrt{10^{\frac{A}{10}}-1}} $$ (which is the relationship between what is normally notated as $\varepsilon$, the variable determining the passband ripple, and $A$ the parameter related to stopband loss), and $$ p_k = 1/\hat{\rho}_k, k=1,2,\ldots,n $$ where $$ \hat{\rho}_k = \sinh(\mu)*real(\rho_k) + j*\cosh(\mu)*imag(\rho_k) $$ where $$ \rho_k = e^{j(\pi(2k-1)/2n + .5\pi)} $$

The poles are given by:

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

 -0.133882765352424 - 0.705787088276200i
 -0.471031461322897 - 0.665351902557852i
 -0.903410457807308 - 0.341931454351974i
 -0.903410457807308 + 0.341931454351974i
 -0.471031461322897 + 0.665351902557852i
 -0.133882765352424 + 0.705787088276200i

And here is the superposition of the magnitude response for each version using Matlab's freqs and zp2tf functions (Matlab in red, Rabiner and Gold in blue---not sure why the phase plot didn't hold as well):

enter image description here

I seem to be missing some normalization factor (k in Matlab-ese, I will refer to it later as $g$) which Rabiner and Gold do not discuss. cheb2ap computes the filter prototype and outputs the zeros, poles and a gain factor given by $$ g = real(\frac{\Pi(-p_k)}{\Pi(-z_k)}) $$ which works out to .01 in cheb2ap and .025 in Rabiner and Gold. However, if I use this in computing the [num,den] = zp2tf(z_rg, p_rg, .025) the magnitude response is still wrong.

enter image description here

It is normalized, but the point at which $|H(\Omega)|^2 = 1/A$ (-40dB) appears to be shifted to the right. Closeup:

enter image description here

  • 2
    $\begingroup$ I changed (and hopefully also corrected) the formula for $\gamma$. I don't have Rabiner and Gold's book, could you maybe add a scan/photo of the relevant page(s)? $\endgroup$
    – Matt L.
    Jul 21 '21 at 10:38
  • $\begingroup$ Ah yes, I messed that one up. It now matches R&B. I don't have access to a scanner at the moment (hurray lockdown!) but I do as you ask when I can. I checked the code and the code computes the correct value for $\gamma$. I also double checked the code vs an implementation in both Matlab and Python so I am 99.9% confident that the discrepancy isn't due to a typo in my C++ code. $\endgroup$
    – dmedine
    Jul 22 '21 at 0:21
  • $\begingroup$ @MattL. Please see my answer below---which I figured out (wait for it) by looking at the figures on the page of text you asked me to upload. $\endgroup$
    – dmedine
    Jul 23 '21 at 2:52

I found out what I was doing wrong. Of course it was my mistake, not the authors of the authoritative text on DSP. I misunderstood the meaning of the parameter $A$ in the Rabiner and Gold formulation. $A$ is not (for example) 40 if you want the cutoff amplitude to be 40dB. Rather, if you want the attenuation in the stop band to be 40dB below unity (or -40dB) then you have this relation: $$ -40 = 20\log_{10}(\frac{1}{A}) $$ in which case $A$ is not $40$, but rather $100$. Using this equation to compute $A$ from the attenuation parameter results in poles that are identical to those produced by the formulae used in Matlab.

  • $\begingroup$ Yes, that makes sense! You can accept your own answer if you like. $\endgroup$
    – Matt L.
    Jul 23 '21 at 6:50

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.