pole/zero-plot --> coefficients?

[Henke]

How do i (manually) calculate the filter coefficients from a pole/zero-plot? Could anyone explain this to a beginner, I really wanna know how it works. (I know complex maths)

[FEV]

Here you will find a rather easy example on how to do it:
http://www.eumus.edu.uy/docentes/jure/d ... #pole-zero
(you have to scroll down a bit)

For more details look at:
http://ccrma.stanford.edu/~jos/filters/

cheers

[Henke]

Great! Thanks a lot!

[HammieTime]

Very old thread, but bump. I have a problem with converting from pz to coefficients. Currently, as I have access to MATLAB student license, I've used this code to generate the coefficients (from 2nd order TF):

%% Note: zeros/poles(1) = angular radius, zeros/poles(2) = normalized angle from 0 to pi
b(1)=0.33;
if (zeros(2) == 0)
b(2) = -(zeros(1)+zeros(1))*0.33;
b(3) = zeros(1)*zeros(1)*0.33;
else
b(2)=-2*zeros(1)*cos(zeros(2))*0.33;
b(3)=zeros(1)*zeros(1)*0.33;
end
% Poles
a(1)=1.0;
if (poles(2) == 0)
a(2) = -(poles(1)+poles(1));
a(3) = poles(1)*poles(1);
else
a(2)=-2*poles(1)*cos(poles(2));
a(3)=poles(1)*poles(1);
end

The thing is, freqz(b,a) and impz(b,a) give reasonable response. That code simply normalises the impulse response so that it doesn't reach above 1. Anyway, the zplane(b,a) gives the correct location of the poles and zeroes (where I put them in polar coordinates). Now, when I use this code in my Synth, everything falls to pieces and it clips like mad. I don't understand, since the filter I've been using now (with the classical bilinear transform method) works just fine, and the freqz, impz and zplane plots are very similar to each other.

My C++ code looks like this, maybe I'm just missing something:

void myBiquad<T>::calcCoeffPolarImpNorm(T* poles, T* zeros)
{
b[0] = 0.33f;
if (zeros[1] == 0) {
b[1] = -(zeros[0]+zeros[0])*0.33f;
b[2] = zeros[0] * zeros[0] * 0.33f;
}
else{
b[1] = -2*zeros[0]*cos(zeros[1])*0.33f;
b[2] = zeros[0]*zeros[0]*0.33f;
}

a[0] = 1.0f;
if (poles[1] == 0) {
a[1] = -(poles[0]+poles[0]);
a[2] = poles[0] * poles[0];
}
else{
a[1] = -2*poles[0] * cos(poles[1]);
a[2] = poles[0] * poles[0];
}
}

I'm using the DF2-transformed method with filter:

T output = *(input+i)*b[0] + xPrev[channel][0] * b[1] + xPrev[channel][1] * b[2] - yPrev[channel][0] * a[1] - yPrev[channel][1] * a[2];

[matt42]

Hi there,

What kind of variable is T in your code?

[HammieTime]

Hi there,

What kind of variable is T in your code?

Hi. It's a template class and now I'm using a float. So the biquad is constructed with myBiquad<float> Biquad1 etc.

[matt42]

OK. Poles and zeros are complex numbers and can't be represented by a float variable.

Besides that for a second order filter if you have the transfer function you should be able to derive the difference equation. There isn't really a need to play with pole zero locations for a synth filter.

[HammieTime]

OK. Poles and zeros are complex numbers and can't be represented by a float variable.

Besides that for a second order filter if you have the transfer function you should be able to derive the difference equation. There isn't really a need to play with pole zero locations for a synth filter.

Neither pole or zero is represented as a complex_t variable since it isn't necessary. For example I can give poles[2]={0.5,pi/4} argument and it simply uses 0.5 as the radius and pi/4 as the angle. This should give a lowpass filter at fs/8 cutoff if the zeros[2]={1.0f,pi}.

[Music Engineer]

i think, the easiest way is to just multiply out the pole/zero product form of the transfer function. in case of a biquad, that should come out as:

Code: Select all

            (1-q1/z) * (1-q2/z)         1 - (q1+q2)/z + q1*q2/z^2
H(z) = k * --------------------- = k * ---------------------------
            (1-p1/z) * (1-p2/z)         1 - (p1+p2)/z + p1*p2/z^2

where k is an overall gain factor (which can taken to be one, if you have none), q1,q2 are your two zeros and p1,p2 are your two poles (this notation is consistent with julius smith's book). from the right hand side, you can just read off the biquad coeffs (b0=k, b1 = -k*(q1+q2), ...). i have no idea, where the 0.33 in your formulas should come from - but if anything, it should probably mean 1./3 - no reason to truncate to two decimal digits in code

[matt42]

Edit:

Sorry, for the constant editing. I think all that needs to be done is normalize the feedforward coefficients so that DC gain = 1.0, at least in the case of a low pass filter. Also, depending on input level/frequency content and resonance you may still get clipping.

In case of 2 zeros at 1, pi it's easy enough to normalise:

Code: Select all

x = (a[0] + a[1] + a[2])/4
b[0] = x
b[1] = 2 * x
b[2] = x

[HammieTime]

Where does the "gain" come from? I mean if the pole radius is close to 1.0f, obviously it's trying to pull the response to infinity at that frequency, but how do you find out the amount of that peak? I think I might have to re-read some of the materials. Intuitively it would mean that if the transfer function hits a "pole", the gain would be 1/(1-pole_radius), if the zero is not close to the pole. Then a radius of 0.9 would lead to a gain of 10 or 20 dB.

The filtering function I've been using now looks like this, maybe not the most efficient but I figured out that by skipping the if-statements, there would be a slight speed increase:

void myBiquad<T>::processChannel(float* channelWrite,const float* channelRead,int channel, int numSamples)
{
for (int i = 0; i < numSamples; i++)
{
T val = *(channelRead+i)*b[0] + xPrev[channel][0] * b[1] + xPrev[channel][1] * b[2] - yPrev[channel][0] * a[1] - yPrev[channel][1] * a[2];
xPrev[channel][1] = xPrev[channel][0];
xPrev[channel][0] = *(channelRead + i);
yPrev[channel][1] = yPrev[channel][0];
yPrev[channel][0] = val;
*(channelWrite + i) = val;
}
}

[HammieTime]

Dual post

[mystran]

Where does the "gain" come from? I mean if the pole radius is close to 1.0f, obviously it's trying to pull the response to infinity at that frequency, but how do you find out the amount of that peak?

You can compute the response at a particular frequency by substituting z=exp(i*w)=cos(w)+i*sin(w) where w=2*pi*f/fs into the transfer function. The complex magnitude gives the gain, the complex angle gives the phase. Note that for DC z=1 and everything stays real (as long as coefficients are real).

[matt42]

I mean if the pole radius is close to 1.0f, obviously it's trying to pull the response to infinity at that frequency, but how do you find out the amount of that peak?

you can find it as mystran explained. The only issue is, depending on the magnitude of the pole, DC gain may be higher than the centre frequency. In those cases normalising the centre frequency will boost frequencies below that point above unity gain

[earlevel]

you can find it as mystran explained. The only issue is, depending on the magnitude of the pole, DC gain may be higher than the centre frequency. In those cases normalising the centre frequency will boost frequencies below that point above unity gain

Yes, and you need to check at the other end too (DC and Nyquist frequency, 0 and pi on the unit circle).

I have a dsp.stackexchange answer here, with way too much detail including code:

https://dsp.stackexchange.com/questions ... 6229#36229