MATLAB's bilinear performs the following steps for a system in zero-pole-gain form

  1. If fp is present, it prewarps:
fp = 2*pi*fp;
fs = fp/tan(fp/fs/2)

otherwise, fs = 2*fs.

  1. It strips any zeros at $±∞$ using
z = z(finite(z));
  1. It transforms the zeros, poles, and gain using
pd = (1+p/fs)./(1-p/fs);    % Do bilinear transformation
zd = (1+z/fs)./(1-z/fs);
kd = real(k*prod(fs-z)./prod(fs-p));
  1. It adds extra zeros at -1 so the resulting system has equivalent numerator and denominator order.

My question is about this line:

kd = real(k*prod(fs-z)./prod(fs-p));

How to derive it? I think it should normalize the gain at $s=0$ and $z=1$, but I don't understand what this line means.

Since it lets fs = 2*fs if fp is not present, it should reads $$ k_d = \mathcal{Re} \{ k_a \frac{\prod_i (2f_s - z_i)}{\prod_i (2f_s - p_i)} \} $$


2 Answers 2


OK, from this answer we know it should be normalized at DC, i.e., let $s=0$ and $z=1$ $$ k_a \frac{\prod_{n=1}^N{(s-z_{a, n})}}{\prod_{n=1}^N{(s-p_{a, n})}}\Bigg|_{s=0} = k_d \frac{\prod_{n=1}^N(z-z_{d, n})}{\prod_{n=1}^N(z-p_{d, n})} \Bigg|_{z=1} $$

$$ k_a \frac{\prod_{n=1}^N{(0-z_{a, n})}}{\prod_{n=1}^N{(0-p_{a, n})}} = k_d \frac{\prod_{n=1}^N(1-z_{d, n})}{\prod_{n=1}^N(1-p_{d, n})} \tag{1} $$

And we know the zero-pole mapping relationship between analog and digital domain $$ z_d = \frac{2f_s+z_a}{2f_s-z_a} \tag{2} $$

$$ p_d = \frac{2f_s+p_a}{2f_s-p_a} \tag{3} $$

Now substituting Eqs. (2) and (3) into (1) yields $$ k_a\prod_{n=1}^N \frac{-z_{a, n}}{-p_{a, n}} = k_d \prod_{n=1}^N \frac{-2z_{a, n} (2f_s-p_{a,n})}{-2p_{a, n}(2f_s-z_{a, n})} $$

and thus $$ k_d = k_a\prod_{n=1}^N \frac{2f_s-z_{a,n}}{2f_s-p_{a, n}} $$


Since z and p are already calculated, it follows that the result of the prod() is a number. Evaluation at z=1 is because z=exp(-i*0)=1, which means the terms are simply added. But, orders higher than 2 are, typically, split into 1st or 2nd order stages. Here it's 1st order, with each stage being made of the respective pole or zero, so what results is the product you see.

In theory, the poles and zeroes are complex conjugate which, when used in a product/sum, the result is a real number. However, given the numerical limitations there will always be residues which will make the result of that product have a non-zero imaginary part, even if it's comparable to machine precision (e.g. 1e-16 or so). The number will still count as complex, so real() is added to cure the infestation.

  • $\begingroup$ Thank you for your answer. I understand the product and the real symbol, but I can't figure out why it's prod(fs-z)./prod(fs-p). I think it's some kind of mathematical derivation. $\endgroup$
    – DSP novice
    Commented Jul 31, 2022 at 10:01
  • $\begingroup$ @DSPnovice Why do you say fs=2*fs? Do you see that in the Matlab function? (I don't have Matlab) It seems strange since, if the user enters fs=2 then 2 Hz should be, not 4. Maybe it's more than meets the eye? $\endgroup$ Commented Jul 31, 2022 at 10:27
  • $\begingroup$ from the first step: If fp is present, it prewarps: fp = 2*pi*fp; fs = fp/tan(fp/fs/2), otherwise, fs = 2*fs. $\endgroup$
    – DSP novice
    Commented Jul 31, 2022 at 11:13
  • $\begingroup$ I think it's only a notation of bilinear constant 2/T for simplicity. $\endgroup$
    – DSP novice
    Commented Jul 31, 2022 at 11:17
  • $\begingroup$ @DSPnovice it doesn't make sense for me. it seems to double the sampling frequency behind your back. If it's an alias for 2/T then why not using it for everything else? Why only there? Well, there are things I don't understand i nthere so please disconsider this answer, I'll delete it after you've seen this reply. $\endgroup$ Commented Jul 31, 2022 at 11:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.