Octave - freqz strange results

I wrote octave code to find transmittance based on zeros and poles:

sps = 44100;     % sample rate [samples/s]

function [b, a] = dbtf(sps, z, p)
b = [ 1 ];
a = [ 1 ];
for i = 1:size(z)(1)
f = z(i, 1);
g = z(i, 2);
b = conv(b, [1, -2*g*cos(2*pi*f/sps), g^2]);
endfor

for i = 1:size(p)(1)
f = p(i, 1);
g = p(i, 2);
a = conv(a, [1, -2*g*cos(2*pi*f/sps), g^2]);
endfor
endfunction

z = [ 0 1; 30 1;            45.8 1;          66 1; ];
p = [           41.7 0.9999;       60 0.9999;     95.2 0.9999; ];
[b, a] = dbtf(sps, z, p);
[h,w] = freqz(b, a, 128000);
semilogx(w/pi+0.00001, abs(h));

It just calculates conjugate zeros and poles of numerator and denominator. Next zeros and poles are added by polynomial multiplying using conv function. Z and p vectors holds frequency/radius pairs.

For such defined zeros and poles I have strange result:

z = [ 0 1; 30 1;            45.8 1;          66 1; ];
p = [           41.7 0.9999;       60 0.9999;     95.2 0.9999; ]; but if I remove zero at 30Hz, the result is as expected:

z = [ 0 1;                  45.8 1;          66 1; ];
p = [           41.7 0.9999;       60 0.9999;     95.2 0.9999; ]; The question: Is it just matter of number precision in octave or my calculating method is wrong ? Alternative question: Could somebody try it on matlab (not octave) and say if results are the same ?

I of course expect notches at frequences of zeros and resonances at frequences of poles. Only second result shows 3 resonances and 3 notches as expected. In the first example frequency response looks completly random (nothing resembles expected notches and resonances).

What you are seeing are indeed numerical artefacts. The conversion from polynomial roots to polynomial coefficients (and vice-versa) is known to be an ill-conditioned problem. In your case the problem occurs at a relatively low polynomial order because zeros and poles are all clustered around the same point (DC), and the poles and zeros almost cancel each other. By removing one zero, as you tried, the problem's numerical condition is improved, simply because the polynomial order becomes smaller. Another way to avoid the problem is to distribute the poles and zeros more evenly around the circle in the complex plane. Note that the available frequencies extend up to $$f_2/2=22050\text{ Hz}.$$ That means that just decreasing the sampling frequency would also improve the numerical condition.

Apart from the numerical effects, your approach for converting zeros and poles to numerator and denominator coefficients is correct. However, it can be done in a much simpler way by using the function zp2tf, which transforms given poles and zeros to transfer function coefficients.

Also note that with your approach you generate a double zero at zero frequency because you produce complex conjugate pairs, even for real-valued zeros.

Finally, the filter you computed is not causal and stable. This has nothing to do with the effects you see, but I believe it's still worth mentioning. So you wouldn't be able to implement that filter with the coefficients you computed. The reason for this is that the transfer function is not proper, i.e., its numerator order is greater than the order of its denominator. This means that apart from the given poles, there occur a number poles at infinity, and that number equals the difference between the numerator and denominator orders. This problem can be solved by adding poles at $$z=0$$, i.e., by multiplying the denominator by $$z^K$$ where $$K$$ is the difference between the polynomial orders of numerator and denominator. This will not change the frequency response, it will just make the filter causal and stable.

EDIT:

For maximum precision you shouldn't compute the coefficients of the numerator and denominator polynomials. Instead you should work with second-order sections. Below is a simple example in Octave/Matlab. Note that I didn't optimize the sequence of sections, neither did I optimize the pole-zero pairing per section. But, as you can see in the figure, the numerical artefacts are greatly reduced.

fs = 44100;
rz = 1; rp = 0.9999;
phi_z = 2*pi*[30; 45.8; 66] / fs;
phi_p = 2*pi*[41.7; 60; 95.2 ] / fs;

m = length(phi_p);
L = 2^15;

a = zeros(m,3); b = zeros(m,3);
H = zeros(m,L);

for k = 1:m,
b(k,:) = [1, -2*rz*cos(phi_z(k)), rz^2];
a(k,:) = [1, -2*rz*cos(phi_p(k)), rp^2];
Htmp = freqz( b(k,:), a(k,:), L );
H(k,:) = Htmp(:).';
end

% add real-valued zero at DC
[H0,w] = freqz( [1, -rz], 1, L );
H = [ H0(:).'; H ];

Htot = prod(H);
semilogx( w/pi+1e-5, abs(Htot) );
axis([1e-5,.01,0,.0005]) • I need resonances at specified frequences (ultimately, I need about 20 resonances). I also can't change sampling frequency. How can I obtain proper polynomial coeficients for that filter? Dec 9 '18 at 22:19
• @apaluch: Please see my edited answer. Dec 10 '18 at 10:31