# Different formulations of Chebyshev II resulting in different poles?

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):

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.

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

• 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)? Jul 21, 2021 at 10:38
• 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. Jul 22, 2021 at 0:21
• @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. Jul 23, 2021 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.