1
$\begingroup$

I'm looking for the practical way to compute the parameters ($a_1$, $a_2$) of a digital filter to match a certain frequency response. I'm studying the 12-poles filter of the vintage component SP0256 which is used in speech synthesis, using bandpass filters to create frequency response of allophones. The 12-poles filter is composed of six 2-poles filters :

$$H(z)=\frac{1}{1-2 a_1 z^{-1}-a_2 z^{-2}}$$

I came accros a github of the FPGA implementation of the SP06256 which provide a matlab example of this filter for the 'EH' allophone. I ploted the 6 frequency responses in blue and the final one in red with matlab :

enter image description here

I could find the formant frequency for the other allophones, but I can't figure out how to compute the filter parameters to match the desired frequency response.

$\endgroup$
  • 1
    $\begingroup$ Is that SP06256 or SP0256? Can you give a link to the FPGA implementation? How is the 12-pole filter composed from the six 2-pole filters? Where are the frequency responses that you want to match to? $\endgroup$ – Olli Niemitalo Dec 9 '18 at 16:38
  • $\begingroup$ github.com/trcwm/Speech256 I cant find the link between filter parameter and the cutoff frequency this is to much computation. I wonder how to determine relation between frequency cutoff and filter parameters. (z=e^jwt and finding where the magnitude = 1/sqrt(2) is the way to go ?) $\endgroup$ – Robert Dawson Dec 9 '18 at 17:16
  • $\begingroup$ There are 2 coefficients per component filter, so a formant frequency is not enough information. You also need something like gain at each formant (peak) frequency. Then again, each component filter will affect the frequencies and gains of the other formants as well, so you may need to optimize the full frequency response, not just one formant at a time. I wonder if there is also some control of the over-all gain of the filter or if the numerator of the transfer function of the cascaded composite filter is really $1.$ $\endgroup$ – Olli Niemitalo Dec 9 '18 at 18:41
  • $\begingroup$ @OlliNiemitalo The numerator can be one with resulting gain, it is just in that case the denominator will be less than one when z is at the frequency of resonance. Each component filter should not effect the other filters; I assume they are all run in parallel so it is just a matter of positioning the pole relative to the unit circle. $\endgroup$ – Dan Boschen Dec 9 '18 at 18:46
  • $\begingroup$ @DanBoschen According to Wikipedia: "The SP0256 realizes its 12-pole filter with a series of cascaded 2-pole IIR filter sections." The frequency response of a 2-pole component filter is not exactly flat at the other formant frequencies than the one it primarily deals with, thus the other formant frequencies will be shifted somewhat. $\endgroup$ – Olli Niemitalo Dec 9 '18 at 18:57
2
$\begingroup$

TL/DR:

For a 2nd order transfer function of a real digital filter given as:

$$H(z) = \frac{1}{1-2A_1z^{-1}-A_2z^{-2}}$$

The filter coefficient $A_1$ and $A_2$ are determined from the gain and resonant frequency to be:

$$A_1 = Kcos(2\pi f_r/f_s) = Kcos(\omega_r)$$

$$A_2 = -(K^2)$$

Where:

$f_r$: Resonant frequency in Hz

$f_s$: Sampling frequency in Hz

$\omega_r$: Normalized angular frequency ($2\pi f_r/f_s$)

$K$: Magnitude of each pole

$G$: Gain, with the gain in dB as $20Log_{10}(G)$

K and G are related with the following formula:

$$G = \frac{1}{(1-K)\sqrt{(1-2K cos(2\omega_r)+K^2)}}$$

I was not able to solve for $K$ as a function of $G$ and $\omega_r$ in closed form (however see reasonable approximation in Update section below); For a given $G$ and $\omega_r$ I would determine $K$ iteratively, or solve for the root of the following equation:

$$(1-K)\sqrt{(1-2K cos(2\omega_r)+K^2)} - \frac{1}{G} = 0$$

Where the solution for K given G is the real positive root between 0 and 1 of

$K^4 + a K^3 + b K^2 + a + c = 0$

with

$a = -2(1+ cos(2\omega_r))$

$b = 4cos(2\omega_r)+2 = -2(1 + a)$

$c = 1- \frac{1}{G^2}$

Thus the solution for K from G can be found from Matlab

roots([1 a b a c])

or Python (NumPy) using:

np.roots([1, a, b, a, c])

Note as Olli pointed out in the comments, the actual resonant frequency will vary slightly from $f_r$; for more on that see his detailed analysis at Precise Centre frequency of an All-pole digital filter

UPDATE

A reasonable approximation for $K$ given $\omega_r$ and $G$ that applies for large gains was determined by @Travis on the mathematics site at this link: https://math.stackexchange.com/questions/3035346/solving-for-ellipse-parameters-given-a-radius-and-angle-challenge-2/3035437#3035437

Resulting in the following for $\omega_r$ not approaching $0, \pi$:

$$K \approx 1 - \frac{1}{Gsin(\omega_r)}$$

With $\omega_r$ = $0, \pi$

$$K = 1- \frac{1}{\sqrt{G}}$$

DETAILS:

The basic 2nd order digital bandpass filters can be constructed as follows:

Position the poles on the z-plane to be close to the unit circle at the desired frequency location for resonance. For a real signal the poles will come in complex conjugate pairs, hence the 2nd order sections. There will also be two trivial zeros at the origin. This is demonstrated in the graphic below showing the pole locations on the z-plane. For the z-plane, the unit circle (drawn) represents the frequency axis, starting with DC at z = 1, and rotating to half the sampling rate at z = -1.

poles on z plane

The two poles shown are located at at 30 degree angle or $2\pi \frac{30}{360}$ in radians. So if the sampling rate was 10 KHz (for example) this would represent a resonant frequency of $10e3 * \frac{30}{360}$ since the 360° around the complete circle represents the sampling rate. The closer the pole is placed to the unit circle, the tighter the resonance and the higher the gain. The gain is determined using the following relationship once doing the math starting with the pole locations in this case of complex conjugage poles:

$$H(z) = \frac{z^2}{(z- p_1)(z-p_2)}$$

$z_1$ and $z_2$ are complex conjugate poles, placed to be close to the unit circle at the angle that represents the resonant frequency:

$p_1 = Ke^{j\omega_r}$

$p_2 = Ke^{-j\omega_r}$

Substituting these in the formula for H(z) along with use of Euler's Identity ($2cos\omega_r = e^{+j\omega_r}+ e^{-j\omega_r}$) results in:

$$H(z) = \frac{1}{1 - 2(Kcos\omega_r) z^{-1} + K^2 z^{-2}}$$

where $K$ is the magnitude of the pole (close to but less than 1), and $\omega_r$ is the angle of the pole in radians, representing the frequency of the resonance as given by the normalized angular frequency.

From this it is clear that:

$$A_1 = Kcos(\omega_r)$$ $$A_2 = -K^2$$

The overall gain for the 2nd order filter is related to the inverse of the distance from the unit circle when z is at resonance to the pole. The exact relationship is complicated given the two pole locations, resulting in a gain that varies depending on the resonant frequency for a given pole magnitude as given by:

$$G = |H(z=e^{j\omega_r}| = \frac{1}{(1-K)|(1-K)cos(\omega_r) + j(1+K)sin(\omega_r)|}$$

Note, in case it helps illuminate a closed form solution for K, that the second part of the denominator has a magnitude that starts at $(1-K)$ for $\omega_r = 0$ and grows as an ellipse to $(1+K)$ at $\omega_r = \pi/2$

The simplest I could get this to is as follows:

$$G = \frac{1}{(1-K)\sqrt{(1-2K cos(2\omega_r)+K^2)}}$$

When the resonant frequency is on the real axis the gain is maximum for a given K:

$$G(\omega_r = 0, \pi) = \frac{1}{(1-K)^2}$$

And when the resonant frequency is on the imaginary axis the gain is at a minimum for a given K:

$$G(\omega_r = \pi/2) = \frac{1}{(1-K)(1+K)}$$

Where K is the magnitude of each pole and G is the resulting gain of the second order filter. (20Log(G) in dB)

The exact solution for K given G for all cases of resonant frequencies is quite entailed, so an iterative solution may be more practical (vary K to converge on G).

For example for the first curve in the OP's plot, from the graph it appears to be resonant at 600 Hz with a gain of 33 dB with an assumed sampling rate of 10KHz. From this we would get:

The resonant normalized angular frequency in radians is:

$$\omega_r = 600/10e3 * 2\pi = 0.377$$

33 dB converted to magnitude:

$$10^{(33/20)} = 44.67 = G$$

The magnitude for each of the two poles is then found from:

$$44.67 = \frac{1}{(1-K)\sqrt{(1-2K cos(2*.377)+K^2)}}$$

$$= \frac{1}{(1-K)\sqrt{(1 - 2K(0.729) + K^2)}}$$

Solving for the root of

$$(1-K)\sqrt{(1 - 1.458 K + K^2)}-1/44.67 = 0$$

Using the approach in the introduction:

$a = -2(1 + cos(2\omega_r)) = -3.4579$ $b = -2(1 + a) = 4.9158$ $c = 1-1/G^2 = 0.9995$

>>> np.roots([1, a, b, a, c])
[0.7294172+0.68505263j 0.7294172-0.68505263j 1.02993616+0.j 0.96914237+0.j]

From which we get from the real root:

$$K = 0.9691$$

Thus $A_1$ and $A_2$ for this case are:

$$A_1 = Kcos(\omega_r) = 0.901$$ $$A_2 = -K^2 = -0.939$$

$$H(z) = \frac{1}{1-2A_1z^{-1}-A_2z^{-2}}$$

= $$\frac{1}{1-1.802z^{-1} + .939z^{-2} }$$

With the frequency response as given below:

Frequency Resposne

$\endgroup$
  • $\begingroup$ That's an excellent analysis @OlliNiemitalo! $\endgroup$ – Dan Boschen Dec 9 '18 at 20:29
  • $\begingroup$ @DanBoschen I wish they had taught me that in engineer school ! Very valuable answer, could you please add more steps between 1st and 2nd equation. First equation do not have any z variable, I don't get it. I'm having trouble understanding identification between A1 and Kcos(phi), A2 and -(K^2) $\endgroup$ – Robert Dawson Dec 9 '18 at 21:53
  • $\begingroup$ @RobertDawson Your welcome. Where are you confused exactly (specifically in the details section that I wrote)? I think I see what may help, adding some detail now... $\endgroup$ – Dan Boschen Dec 9 '18 at 22:01
  • $\begingroup$ @DanBoschen At H(z) = 1/((1-z1)(1-z2)), shouldn't there be the z variable, something like H(z) = 1/((z-p1)(z-p2)) p1,p2 representing the 2 poles ? $\endgroup$ – Robert Dawson Dec 9 '18 at 22:05
  • $\begingroup$ @RobertDawson ah yes, of course-- sorry about that confusion. I fixed it, let me know if it now makes sense. (I also added the $z^2$ that is supposed to be in the numerator, which then comes out to your equation given in negative exponents of z; dividing numerator and denominator by $z^2$ will give the equivalent result. $\endgroup$ – Dan Boschen Dec 9 '18 at 22:08

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.