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

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.

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

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.

• 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. Jul 31, 2022 at 10:01
• @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? Jul 31, 2022 at 10:27
• from the first step: If fp is present, it prewarps: fp = 2*pi*fp; fs = fp/tan(fp/fs/2), otherwise, fs = 2*fs. Jul 31, 2022 at 11:13
• I think it's only a notation of bilinear constant 2/T for simplicity. Jul 31, 2022 at 11:17
• @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. Jul 31, 2022 at 11:24