Bilinear transform for a system in zero-pole-gain form

Question

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.

2. It strips any zeros at $±∞$ using

z = z(finite(z));

3. 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));

4. 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)} \} $$

Answer 1

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}} $$

Answer 2

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.