# Bilinear transform with pre-warping for systems other than classical filters

I only seem to be able to find online information about applying bilinear transform + pre-warping to filters (like butterworth, etc.) with only one edge frequency that is purposely 'designed' into it. So the BT is easy to apply, just substitute with the usual tangent formula including edge frequency and sampling frequency.

But what if my analog rational transfer function describes a physical system (like an elastic mechanical body with several eigen-frequencies, or an electric circuit network)? Do I have to apply the BT to every pole/zero separately then? And what if the transfer function comes from measurements - I don't know the poles/zeros, because I did not 'design' the 'filter'? Does it really mean I have to factorize it (which we know is sensitive to error) just to transform it from analog to digital?

PS, Sidenote: I am going to implement the numerics involved (in C#) myself. So Matlab is certainly a nice thing, but it does not help me if not all relevant code is publicly available.

• that "factorize it" is, i presume, breaking into partial fractions for the purpose of using the method called Impulse Invariant in lieu of Bilinear Transform. For me, the question regarding this choice has always been: "Which domain is most important to you?" If you wanna match time-domain behavior, use Impulse Invariant. If you wanna match frequency-domain behavior, use Bilinear Transform and pre-warp every design frequency that you have independent control over. If you have a high-order filter, there will likely be more than one frequency that you can pre-warp to match. May 17, 2019 at 23:40
• @robertbristow-johnson: my feeling is that the impulse invariant method is indeed favorable in my case. The problem I have with understanding this method is: partial fraction decomposition requires finding zeros and poles numerically. AFAIK general root finding is pretty tough in itself, and it gets even tougher if there are (presumably) removable singularities in the transfer function, because of roundoff errors. How can I avoid ending up with a mess of code with high error sensitivity when I try to implement impulse invariance myself? Is there an open-source code base for this? May 19, 2019 at 18:58
• Octave is a free Matlab clone, so you could check if the functions you need exist and use their code to understand the methods. I know that the function invfreqz.m exists in the Octave signal processing toolbox. May 19, 2019 at 19:13
• Matt, thanks for the hint. I actually have already tried to do so with the example of the c2d-function you mentioned (c2d.m). Apparently I seem to fail to understand what this m-file means. There is some obscure pseudo-code in it, each line starting with "%!". Nothing I am familiar with from Matlab. May 19, 2019 at 19:29
• @oliver: These are just test blocks, nothing to do with the function code itself. May 19, 2019 at 20:11

The bilinear transform ("Tustin's method") is indeed mainly used for transforming frequency selective filters with magnitude responses that are optimal with respect to some criterium, such as Butterworth, Chebyshev, or Cauer filters.

For more general systems, the bilinear transform is usually not the best choice because of the frequency warping, which mainly affects frequency domain behavior close to Nyquist. For frequency-selective filters, cut-off frequencies can be pre-warped such that the cut-off frequencies of the resulting discrete-time system are correct, but for more general systems it is usually not clear which frequencies should be matched and which shouldn't. Also, for filters with approximately piece-wise constant magnitude responses, frequency warping is not a big problem as long as the cut-off frequencies are mapped correctly. For general systems one would usually prefer the error between the continuous-time system and the discrete-time system to be distributed more evenly in the frequency domain, which is not possible with the bilinear transform.

There are several methods for converting continuous-time systems to discrete-time systems. All of them have advantages and disadvantages, and it depends on the application which one is to be preferred. Take a look at Matlab's c2d function, which implements several methods. Also take a look at this question and its answers.

For approximating a measured response, you would directly design the corresponding system in the discrete-time domain. For general frequency responses, the frequency sampling method is very simple and effective. If an IIR approximation is desired, the Matlab function invfreqz could be used. It tries to fit a rational transfer function to given frequency response data.

• I guess the frequency sampling method is out of the picture when it comes to physical systems, because they more often than not have infinite impulse response (vibration, exponential falloff). May 18, 2019 at 20:35
• @oliver: I wouldn't say so, because you can approximate any (stable) system by a system with a finite impulse response. May 19, 2019 at 8:51
• Matt, I think the important thing is, how approximation is defined in this case. The fact that an excited mechanical system comes to rest in finite time might be practically irrelevant if that duration is long enough, because any measurement is finite as well. So it is mainly about whether I want to extrapolate the system evolution beyond the measurement. But if I have a FIR with say 10000 samples, it gets a little impractical to perform and I will be better off transforming the continuous transfer function to fourier space in the first place. May 19, 2019 at 10:44
• @oliver: It depends on a lot of things that are not stated in your question: memory of the system, sample rate, etc.; FIR is an option in general, it might not be in your case, but we can't know. And what is practical and what not depends on the available computational power, also something we can't know. May 19, 2019 at 11:33
• @oliver: Your question is OK, if you don't mind getting more general answers due to lack of specific information. I've added a few sentences to my answer referring to IIR approximation of given frequency response data. May 19, 2019 at 13:33

If you want to convert a physical system transfer function the most "accurate" way is usually to convert it by use the "step-invariance" method where you add the effect of your zero-order hold to the Laplace transform of your process and the you convert it to the Z domain. It correctly models the zero-order hold effect of your DAC.

A second method that works well is to design your controller in the s-domain and to convert it in the z-domain using the bilinear transform. It works well when your sampling rate is 20 times fast or more than the open loop bandwidth. You can also consider modeling your zero-order hold (and all the other delays) using the Padé method for more accuracy.

• As to the first paragraph, I have to admit that I don't know remotely what you mean. Sorry for that, I'm a physicist with only a modest signal processing background. As to the second paragraph: well yes, I knew that pre-warping is mainly for features close to Nyquist, so if I am far from it, I can refrain from it. My question was intended for systems that might be considered being on the fringe. May 18, 2019 at 20:21
• I don't think you should focus on the pre warping. I didn't mention it in my second paragraph
– Ben
May 18, 2019 at 20:51
• Sure you didn't mention it, but pre-warping (as it contains "edge" frequency) is exactly what makes it difficult for me to translate the most general analog system (like a damped harmonic oscillator with a certain resonance frequency and width aka damping) more accurately to a digital system. If I am willing to accept deviation of resonance frequency near Nyquist, then, of course I am fine with the simple bilinear transform. But at the moment I am still surveying what I can find about step-invariance. May 18, 2019 at 20:56
• What is your sampling frequency and what is the frequency of your signal of interest?
– Ben
May 18, 2019 at 21:18
• That may depend... ;-) What I am trying to say is: it could be a relatively poorly sampled signal and then the resonance may be pretty close to Nyquist. Or it may be a well sampled signal where the resonance is comfortably far away from Nyquist. Of course I could always upsample, but I'd prefer a one-stage process for going from analog to digital. May 18, 2019 at 21:23