1
$\begingroup$

I'm not an engineer and have essentially taught myself all I know - this present problem is giving me some problems.

I have an all-pole filter (a gammatone filter - impulse response is essentially a damped cosine) implemented as a cascade of identical 2nd-order pole pair sections. The poles are located in the z-domain at $\exp(-\beta + i\theta)$: $\beta$ is the bandwidth in radians/sample and $\theta$ is the 'ringing' frequency in radians/sample.

Given a value for $\theta$, this pole pair does not produce a filter with peak gain at the frequency $\theta$ - I'm trying to find out exactly what the offset is (i.e. at what frequency a filter with a given $\theta$ actually has peak magnitude), so that I can produce filters that peak at precisely a given arbitrary frequency.

I've found an equation for calculating what this offset is given poles located in the s-domain at $-\beta + i\theta$:

$$\mathrm{f_{c,actual}} = \theta \sqrt{1 - \dfrac{\beta^2}{\theta^2}}$$

This formula provides results consistent with observations up to a certain frequency, though I believe it needs to be modified to produce results that work for arbitrary discrete spectrums.. I've tried a number of things and just can't quite figure it out.

Any help is greatly appreciated.

$\endgroup$
8
$\begingroup$

For a filter consisting of a complex conjugate pair of poles, the $z$-domain transfer function is: $$H(z) = \frac{a}{\left(1-r(\cos{\theta}-i\sin{\theta})z^{-1}\right)\left(1-r(\cos{\theta}+i\sin{\theta})z^{-1}\right)}\\ = \frac{a}{1 - 2r\cos(\theta)z^{-1} + r^2z^{-2}},$$

where $a$ is a constant by which the magnitude frequency response can be normalized, $r$ is the pole radius, and $\theta$ is the pole angle. The peak of the magnitude frequency response will not be exactly at frequency $\theta$, but at a nearby frequency $\omega_0$ at which the derivative

$$\frac{|a|2r\sin(\omega)\left(2r\cos\omega - (r^2 + 1)\cos\theta\right)}{\left(4r^2\cos^2\theta - 4r(r^2 + 1)\cos(\theta)\cos(\omega) + 4r^2\cos^2\omega + r^4 - 2r^2 + 1\right)^{3/2}}$$

of the magnitude frequency response $|H(e^{i\omega})|$ with respect to $\omega$ equals 0. That frequency is:

$$ω_0 = \operatorname{acos}\left(\frac{(r^2 + 1)\cos\theta}{2r}\right)\text{ if }-1 < (r^2 + 1)\cos(\theta)/(2r) < 1$$

No-go droplets
Figure 1. The condition for that the magnitude frequency response peaks elsewhere than at $w_0=0$ or $w_0=\pi$ is that poles stay outside these droplets.

In the following it is assumed that the condition is satisfied as otherwise the peak would be trivially at frequency $0$ (if $\theta < \frac{\pi}{2}$) or $\pi$, whichever is closer to the pole. The normalized magnitude frequency response with $a$ chosen so that $|H(e^{i\omega_0})| = 1$ reads:

$$|H(e^{i\omega})| = \sqrt{\frac{(r^4 - 2r^2 + 1)\sin^2\theta}{4r^2\cos^2\theta - 4r(r^2 + 1)\cos(\theta)\cos(\omega) + 4r^2\cos^2\omega + r^4 - 2r^2 + 1}}.$$

Given $\omega_0$ and $r$:

$$\theta = \operatorname{acos}\left(\frac{2r\cos\omega_0}{r^2 + 1}\right)\tag{1}$$

which is something you can use directly to keep the peak of your all-pole bandpass filter where you want it:

Constant <span class=$\omega_0$ pole loci">
Figure 2. Pole loci for a selection of constant $\omega_0$ traced over $0 < r < 1$.

Applying the substitution given by equation (1) gives:

$$H(z) = \frac{a}{1 - \frac{4r^2\cos(\omega_0)}{r^2 + 1}z^{-1} + r^2z^{-2}}$$

Just like $\theta$ does not exactly give the peak frequency, $r$ does not give the exact bandwidth. The exact -3 dB to -3 dB bandwidth $B_\text{-3 dB}$:

Exact bandwidth
Figure 3. -3 dB bandwidth of the filter, with $\omega_0 = \frac{\pi}{4}$ and $r=0.9$.

is given by:

$$B_\text{-3 dB} = \operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + (r - 1/r)\sin\theta}{2}\right)\\- \operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + (1/r - r)\sin\theta}{2}\right)$$

This was found by solving the values of $\omega$ for which $|H(e^{i\omega})| = \sqrt{2}/2$. The first $\operatorname{acos}$ gives the solution found to the right from the peak and the latter $\operatorname{acos}$ the one to the left. More generally, for a bandwidth from magnitude frequency response value $g$ to $g$:

$$B_g = \operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + \sqrt{1/g^2 - 1}(r - 1/r)\sin\theta}{2}\right)\\-\operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + \sqrt{1/g^2 - 1}(1/r - r)\sin\theta}{2}\right)$$

This is needed when designing a composite filter, for example a gammatone filter, that is a cascade of identical two-pole filters. Usually we want to design the filter to have a given bandwidth, but it is very difficult to solve $r$ from the above. So we resort to approximation. Using the value $D_2$ of the second derivative of the magnitude frequency response at $\omega_0$:

$$D_2 = \frac{(r^2 + 1)^2\cos^2\theta - 4r^2}{(r^2 - 1)^2(1 - \cos^2\theta)},\tag{2}$$

we can approximate the magnitude frequency response with that of a complex one-pole filter that has at $\omega_0$ the magnitude frequency response and its first two derivatives matched to those of the two-pole filter:

One-pole approximation Figure 4. One-pole approximation with pole angle $\omega_0$.

The magnitude frequency response of the one-pole filter is:

$$|\hat H(e^{i\omega})| = \frac{\sqrt{2}\left(\sqrt{1 - 4D_2} - 1\right)}{2\sqrt{\left(2D_2\cos(\omega - \omega_0) - 2D_2 + 1\right)\left(1 - \sqrt{1 - 4D_2} - 2D_2\right)}}.$$

It is symmetrical around $\omega_0$ and has $g$ to $g$ bandwidth:

$$\hat B_g = 2\operatorname{acos}\left(\frac{1/g^2 - 1}{2D_2} + 1\right),$$

or inversely and more usefully:

$$D_2 = \frac{1 - g}{2g^2\left(\cos(\hat B_g/2) - 1\right)}.$$

Equations $(1)$ and $(2)$ when combined give:

$$D_2 = \frac{2r^2(r^2 + 1)^2\left(1 - \cos(2\omega_0)\right)}{(r^2 - 1)^2(2r^2\cos(2\omega_0) - r^4 - 1)}.$$

Solving that for $r$ gives a cumbersome equation with $\omega_0$ and $D_2$ as arguments that can be used to set the approximate bandwidth of the two-pole filter (c language compatible plaintext with w_0 = $\omega_0$ and d_2 = $D_2$):

r=sqrt(-sqrt(2)*sqrt(sin(w_0))*sqrt((-sin(w_0)*(pow(d_2+1,2)*pow(cos(w_0),2)
+d_2-1)*sqrt(-pow(d_2+1,2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2+1)+pow(d_2+1,3)
*pow(cos(w_0),4)-(d_2+1)*(pow(d_2,2)-d_2+2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2
+1)/sqrt(-pow(d_2+1,2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2+1))+sin(w_0)*sqrt(-
pow(d_2+1,2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2+1)-(d_2+1)*pow(cos(w_0),2)+1)
/sqrt(-d_2)

or with alternative symbols:

r = √(- √2·√(SIN(ω_0))·√((- SIN(ω_0)·((d_2 + 1)^2·COS(ω_0)^2 + d_2 - 1)·√(-
(d_2 + 1)^2·COS(ω_0)^2 + d_2^2 - 2·d_2 + 1) + (d_2 + 1)^3·COS(ω_0)^4 - (d_2
+ 1)·(d_2^2 - d_2 + 2)·COS(ω_0)^2 + d_2^2 - 2·d_2 + 1)/√(- (d_2 + 1)^2
·COS(ω_0)^2 + d_2^2 - 2·d_2 + 1)) + SIN(ω_0)·√(- (d_2 + 1)^2·COS(ω_0)^2 +
d_2^2 - 2·d_2 + 1) - (d_2 + 1)·COS(ω_0)^2 + 1)/√(-d_2)

When $\hat B_\text{-3 dB} = \hat B_{\sqrt{2}/2}$ is set proportional to $\omega_0$ by a constant $c$ stepped from 0.1 to 1.0, the upper half-plane pole resides on tracks like these:

Pole loci for bandwidth proportional to peak frequency
Figure 5. Pole loci for bandwidth proportional to peak frequency with proportionality constant $c$.

$\endgroup$

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.