Im confused by Matlab's dontinuous-to-discrete conversion methods and how they work

Hello community

Recently, I have been working a lot with the transformation from the s-domain to the z-somain using Matlab's c2d function. Currently I'm trying to get an overview about the different discretization methods, i.e. how do they work and when to use them. I'm confused by the notation Matlab uses compared to other literature. Everything I am asking is related to this documentation.

All I am interested in are transfer functions in the continuous-time and discrete-time domain. Meaning that given a transfer function , which steps do I have to perform to get the transfer function . This is what I came up with so far:

Can someone help me complete this table?

Furthermore, there are some more questions associated to this topic

1) As far as I understood, not every transformation is a simple substitution of the variable s. However, is it possible to give a general relationship between the transfer function in the s-domain and z-domain for all cases?

2) According to the documentation, zero-order hold and first-order hold are used best for good matching in the time dimain for "staircase" and "piecewise linear" inputs. I do not know what to do with this information, i.e. why would an analog system be driven by such signals. What is the connection here?

3) Zero-Pole Matching should be used for good matching in the frequency domain according to Matlab's documentation. However, other sources state that this is not the case. I link the related Wikipedia article here (more sources stating this are cited there). In my case I am interested in modelling analog audio components in the digital domain. Hence, precise matching of the frequency response is important. Now I am not sure whether Zero-Pole matching is a good method or not.

4) The Tustin approximation also yields good matching in the frequency domain according to the documentation. From a theoretical point of view this cannot be true in general since the whole frequency axis is wrapped around the unit circle in a non-linear matter, i.e. there must be frequency distortion. The documentation further introduces the concept of the prewrapper frequency, which allows to capture the dynamics of an analog system on a particular frequency. This all yields to the same question as above: If I want to copy the behaviour of an system with audio frequencies (0Hz - 20 kHz), is this method suitable or not?

5) One transformation which I could not find in Matlab is the backwards difference. This transformation is used e.g. in Simulink for the discrete-time derivative block (assuming that a derivative in the continuous domain has the ideal transfer function ). I'm not sure whether this substitution even es a real transformation since I was unable to find related literature. Somebody has more information about it?

6) Why can I not use some transformations for non-causal systems? For example, I want to converte an (ideal) continuous-time derivative:

s = tf('s');
Ts = 5;         % Sampling period
Gs = s;         % Continuous-time TF
Gz = c2d(Gs,Ts,'impulse');  % Discrete-time TF

This gives an error for any method except Tustin's approximation and Zero-Pole matching. I know that there are stability issues, i.e. a contunuous-time derivative is unstable per definition and some transformation map a stable analog system to an unstable discrete system. Is this the problem? And why can Tustin's methor or Zero-Pole mapping handle this?

Answer 1

[Paul]

I'll give an initial answer and then come back and edit if you have any more questions.

As far as the table goes:

The Zero-order Hold method is called just that, i.e., the zero-order hold approximation. It is not the Forward Difference or Euler's method (there was actuallly another post a few weeks ago where this same subject came up). For example:

>> H=tf([1 1],conv([1 2],[1 3]))  % simple transfer function

H =
 
      s + 1
  -------------
  s^2 + 5 s + 6
 
Continuous-time transfer function.

>> Ts = 0.1;
>> Hd = c2d(H,Ts,'zoh')  % zoh approximation

Hd =
 
   0.08215 z - 0.07432
  ---------------------
  z^2 - 1.56 z + 0.6065
 
Sample time: 0.1 seconds
Discrete-time transfer function.

>> s = tf([1 -1],Ts,Ts) % forward difference substitution

s =
 
  z - 1
  -----
   0.1
 
Sample time: 0.1 seconds
Discrete-time transfer function.

>> minreal((s + 1)/(s*s + 5*s + 6))  % sub into original transfer function, not the same as zoh approximation

ans =
 
     0.1 z - 0.09
  ------------------
  z^2 - 1.5 z + 0.56
 
Sample time: 0.1 seconds
Discrete-time transfer function.

I think that the transformation in your third column is correct for the Forward Difference method. As far as I know, the CST doesn't support the Forward Difference or Backward Difference methods, at least not officially (see below) which is suprising to me.

The zoh and foh and impulse invariant and zero-pole matching methods can't be described in terms of a simple substitution between s and z, like can be done for Tustins method. But the zoh and impulse invarant and zero-pole matching methods are pretty straightforward. I'm not as familiar with the foh method.

Impulse invariant mapping is often called impulse invariance.

Reference [2] in the doc page you linked is a very good resource.

As to your other questions:

1. I don't understand what you mean by a "general relationship." Can you clarify?

2. Suppose you have an analog device that is being driven a compute that develops inputs to the device in discrete time. In order to link the device to the computer, you have to convert the discrete time signal into an an analog signal, i.e., a D-to-A converter. A D-to-A converter that maintains a constant signal until the next discrete time update by the the computer is a zero-order-hold. A first-order-hold extrapolates linearly until the next computer update. Hence, the zoh and foh methods, IMO, are primarily appropriate when you need a discrete model of an analog system with a D-to-A at its input. Others may disagree.

3. Zero-pole matching or any of the other non-hold methods may or may not be what you need. You should try them all for your application, plot up responses of interest, and see which result gives you the best answer for your use case.

4. All of the methods will result in some distortion as the frequency approaches the Nyquist frequency, some more than others.The prewarp frequency with Tustin's method yields the perfect match at the prewarp frequency, not over a range of frequencies. An example of when you might want to use pre-warping would be if you had an analog notch filter and you want its discrete apprximation to be perfectly matched at the notch frequency.

5. c2d does not support forward or backward difference, Well, at least not officially. This actually works, though undocumented

c2d(H,T,'forward')

Edit: Actually, it doesn't work. See comment below.

6. I think the short answer here, and I don't mean to be flip, is that the math doesn't work out for those methods. For example, the impulse invariant method is defined in terms of samples of the continuous time impulse response of the continuous time transfer function. The impulse response of H(s) = s is the derivative of the Dirac delta function. How can that function be sampled at t=0?

[Dominik Hiltbrunner]

Hello Paul

Thank you very much for your respons, much appreciated.

First of all I corrected my table which turend out to be wrong.

1. I don't understand what you mean by a "general relationship." Can you clarify?

By that I mean what is the connection between the transfer function in the s-domain and z-domain for each transformation. I expected that there was a closed formula for each transformation that I can use to get the corresponding z-TF, given the s-TF. Apparently, this turned out to be wrong. You already mentioned this avobe:

"The zoh and foh and impulse invariant and zero-pole matching methods can't be described in terms of a simple substitution between s and z, like can be done for Tustins method."

So that quesstion is answered.

The Zero-order Hold method is called just that, i.e., the zero-order hold approximation. It is not the Forward Difference or Euler's method

You're right. The reason why I believed this was true is this thread. However, as it turned out the integrator was just a special case where this substitution yields the same result as the zoh transformation. Trying out some examples in Matlab confirms this.

I think that the transformation in your third column is correct for the Forward Difference method.

Yes, that's what my text book says too (no actual book, just some slides from my university with no source...)

You should try them all for your application, plot up responses of interest, and see which result gives you the best answer for your use case.

From a lot of trial & error I know that the Tustin method works best for me. However, I would like to have some theoretical background to prove that this method is best rather than just a heuristic approach. That's why I try to get an overview about certain transformations.

4. All of the methods will result in some distortion as the frequency approaches the Nyquist frequency, some more than others

My "text book" tells me that e.g. the impulse invariant method mapps the frequency axis in a linear matter around the unit circle. Hence, no distortion. But this is only true for frequencies below the nyquist frequency since this transformation wrapps the frequency axis periodically around the unit circle, hence aliasing. This is also what wiki says.

This actually works, though undocumented: c2d(H,T,'forward')

Very interesting. Sadly,

edit c2d

does not give a deeper insight.

6. I think the short answer here, and I don't mean to be flip, is that the math doesn't work out for those methods.

As you already mentioned, this makes sense e.g. for the impulse invariant transformation whose goal is to match the impulse response of the analog filter with the discrete one. Systems like a simple derivative cannot be described in this way. On the other hand,

c2d(s,0.1,'forward')

throws an error too, although this should work easily. Note that this is note the best example since it is not officially supported.

I will probably also edit my original post and add some information once I got them.

[Paul]

1. The forward, backward, and Tustin transformations can all be derived by different approaches to numerical integration of a continuous signal. A deeper connection lies in the fundamental relationhip z = exp(sT). The forward rule, z = 1 + Ts, is simply the first two terms in the power series expansion of exp(sT). The backward and Tustin rules are two rational function approximations to exp(sT). See this link for more discussion about these approximations. In Matlab, you can get approximations to exp(-sT) (note the negative sign). See

doc pade

that you can use to generate higher order substitutions between z^-1 and s, which you can solve to get s in terms of z and then substitute as usual (I've never actually tried this, and so so don't know if there are any pitfalls waiting to be found, but I don't see why it shouldn't work, at least in theory).

The zoh and impulse and matched (and I suspect foh as well) all also rely on the relationship z = exp(sT). For our example, the impulse method can be replicated as follows:

>> H=tf([1 1],conv([1 2],[1 3]))

H =
 
      s + 1
  -------------
  s^2 + 5 s + 6
 
Continuous-time transfer function.

>> [r,p,k]=residue([1 1],conv([1 2],[1 3]));
>> T=0.1;
>> % note the exp(sT) terms evaluated at the poles of H(s)
>> % this would be more complicated with repeated poles
>> % fixed typo from original post
>> Hd=(tf(r(1)*[1 0],[1 -exp(p(1)*T)],T) + tf(r(2)*[1 0],[1 -exp(p(2)*T)],T))*T

Hd =
 
   0.1 z^2 - 0.08966 z
  ---------------------
  z^2 - 1.56 z + 0.6065
 
Sample time: 0.1 seconds
Discrete-time transfer function.
>> c2d(H,T,'impulse')

ans =
 
   0.1 z^2 - 0.08966 z
  ---------------------
  z^2 - 1.56 z + 0.6065
 
Sample time: 0.1 seconds
Discrete-time transfer function.

Whether or not Tustin's method is the best is really up to you. After all, "best" is subjective.

4. Maybe I don't understand what your text book means by distortion.What I mean is that the discrete approximation starts to diverge from the continuous time transfer function as frequency gets closer to the Nyquist frequency. For example:

>> bode(H,c2d(H,T,'impulse'))

I may have been too hasty in pointing out that forward is usable but undocumented. From a simple example:

>> H=tf(1,[1 1])

H =
 
    1
  -----
  s + 1
 
Continuous-time transfer function.

>> c2d(H,T,'forward')

ans =
 
  0.04837 z + 0.04679
  -------------------
      z - 0.9048
 
Sample time: 0.1 seconds
Discrete-time transfer function.

That result doesn't look correct to me; you can check it yourself using the substitution in your table. Hmm, an undocumented feature that returns an incorrect(?) result.

The reason this result is incorrect is because it's not really the forward rule. As it turns out, it appears that c2d only cares about the first character of the 'method' argument. So specifying 'forward' actually means 'foh'! Examples:

>> isequal(c2d(H,T,'fzdf'),c2d(H,T,'foh'))

ans =

  logical

   1

>> isequal(c2d(H,T,'zebra'),c2d(H,T,'zoh'))

ans =

  logical

   1

6. You're right. I don't think using the forward method is a good example, because it's undocumented and might not even be correctly implemented. I really don't know why c2d doesn't officially support forward and backward. If nothing else, they are good to use for textbook examples.

[Dominik Hiltbrunner]

Ok we're making some progess, which is nice.

First of all let's update our table.

Of course we do not just take this for granted. Let's compare it with Matlabs methods (no guarantee for corectness). I do not verify the simple substitution methods since those are trivial.

Impulse-invariant mapping

% Case 1: Impulse-invarnat transformation
clc; clear; close all;
z = tf('z');
Ts = 0.1;

Hs = tf([1 -1],conv([1 1],[1 2])) % z1 = 1, p1 = -1, p2 = -2;
Hz = c2d(Hs,Ts,'impulse')

[R,P,K] = residue([1 -1],conv([1 1],[1 2]));
N = length(P);

% This code only works in simple cases where
% Gs only has poles with multiplicity 1
Hz_copy = 0;
for n=1:N
    Hz_copy = Hz_copy + R(n)*z/(z-exp(Ts*P(n)));
end
Hz_copy = Hz_copy * Ts

Hz =
 
    0.1 z^2 - 0.1077 z
  ----------------------
  z^2 - 1.724 z + 0.7408
 
Sample time: 0.1 seconds
Discrete-time transfer function.


Hz_copy =
 
    0.1 z^2 - 0.1077 z
  ----------------------
  z^2 - 1.724 z + 0.7408
 
Sample time: unspecified
Discrete-time transfer function.

Pole-Zero mapping

% Case 2: Pole-Zero Mapping
clc; clear; close all;
z = tf('z');
Ts = 0.1;

Hs = tf([1 -1],conv([1 1],[1 2])) % z1 = 1, p1 = -1, p2 = -2;
Hz = c2d(Hs,Ts,'matched')

% Here we need the zeros too, so we use pzmap
[P,Z] = pzmap(Hs);
M = length(Z);
N = length(P);

% Construct numerator and denominator
Num = 1;
Den = 1;
for n=1:M
    Num = Num * (z-exp(Z(n)*Ts));
end
for n=1:N
    Den = Den * (z-exp(P(n)*Ts));
end

% Match DC gain
Hz_copy = Num/Den;
Hs_DC = dcgain(Hs);
Hz_DC = dcgain(Hz_copy);
Hz_copy = minreal(Hz_copy * Hs_DC/Hz_DC)

Hz =
 
   0.08201 z - 0.09063
  ----------------------
  z^2 - 1.724 z + 0.7408
 
Sample time: 0.1 seconds
Discrete-time transfer function.


Hz_copy =
 
   0.08201 z - 0.09063
  ----------------------
  z^2 - 1.724 z + 0.7408
 
Sample time: unspecified
Discrete-time transfer function.

Zero-order hold

% Case 3: Zero-order hold
syms s z n t
Ts = 0.1;

Hs = tf([1 -1],conv([1 1],[1 2])) % z1 = 1, p1 = -1, p2 = -2;
Hz = c2d(Hs,Ts,'zoh')

Hs_sym = poly2sym([1 -1],s)/poly2sym(conv([1 1],[1 2]),s);
h = ilaplace(Hs_sym/s);
h = subs(h,t,n*Ts);
Hz_sym = collect((1-z^-1)*ztrans(h));
[Num,Den] = numden(Hz_sym);
Nc = coeffs(Num,z,'all');
Dc = coeffs(Den,z,'all');
Hz_copy = minreal(tf(double(Nc),double(Dc),Ts))

Hz =
 
    0.08158 z - 0.0902
  ----------------------
  z^2 - 1.724 z + 0.7408
 
Sample time: 0.1 seconds
Discrete-time transfer function.


Hz_copy =
 
    0.08158 z - 0.0902
  ----------------------
  z^2 - 1.724 z + 0.7408
 
Sample time: 0.1 seconds
Discrete-time transfer function.

What is missing is the first-order hold and the least-squares method. So far I do not know how to write down their transformation in one line. I will add them when I have the information.

[Paul]

Forward/Backward - Your table should say "Not Supported"

Impulse Invariance: You should add a comment to your table that the formula you show is only valid for the form of H(s) in the table, i.e., H(s) does not have repeated poles. Repeated poles have to be handled differently.

ZOH: It looks like you got the right answer for this case using the approach with the Symbolic Toolbox. You can also do the computation directly using Control System Toolbox functions:

>> H=tf([1 1],conv([1 2],[1 3]))

H =
 
      s + 1
  -------------
  s^2 + 5 s + 6
 
Continuous-time transfer function.

>> Hd=zpk(tf([1 -1],[1 0],T)*c2d(H*tf(1,[1 0]),T,'impulse')/T)

Hd =
 
   0.082153 z (z-1) (z-0.9047)
  -----------------------------
  z (z-1) (z-0.8187) (z-0.7408)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

>> minreal(Hd)

ans =
 
   0.082153 (z-0.9047)
  ---------------------
  (z-0.8187) (z-0.7408)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

>> zpk(c2d(H,T,'zoh'))

ans =
 
   0.082153 (z-0.9047)
  ---------------------
  (z-0.8187) (z-0.7408)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

Pole-Zero Matching: I don't think your formula is quite right. As you show, N is the order of the demonimator and M is the order of the numerator. After you map the poles and zeros through z = exp(sT), you also need to include an additonal N-M-1 zeros at z = -1 (some people might include N-M zeros at z = -1). The reason is that your original transfer function has N-M zeros at s = j*inf, so your discrete time transfer function should have the same number of zeros at the highest frequency point on the unit circle, which is z = -1. However, there may be a good reason to only have N-M-1 zeros at z = -1, which is what c2d does; consult a good text book for further discussion. Here's an example:

>> H = zpk(-1,-[2 3 4 5],1)

H =
 
           (s+1)
  -----------------------
  (s+2) (s+3) (s+4) (s+5)
 
Continuous-time zero/pole/gain model.

>> Hd = c2d(H,T,'matched')  % N=4, M=1, note that 4-1-1 = 2 zeros are included at z = -1

Hd =
 
         0.00013342 (z-0.9048) (z+1)^2
  -------------------------------------------
  (z-0.8187) (z-0.7408) (z-0.6703) (z-0.6065)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

First Order Hold: The text book approach (at least what I'm familiar with) assumes a causal FOH on the inputs. The transform formula for your table would be:

Hd = (1/T) * (1 - z^-1)^2 * Z(Linv(H(s) * (Ts + 1) / s^2))

This can be implemented directly as:

>> Hd = zpk(1,0,1,T)^2*c2d(zpk(H)*zpk(tf([T 1],[1 0 0])),T,'impulse')/T^2

Hd =
 
  0.12603 z (z-1)^2 (z-0.9049) (z-0.347)
  --------------------------------------
    z^2 (z-1)^2 (z-0.8187) (z-0.7408)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

>> Hd = minreal(Hd)

Hd =
 
  0.12603 (z-0.9049) (z-0.347)
  ----------------------------
    z (z-0.8187) (z-0.7408)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

However, the CST uses a different approach and assumes a non-causal FOH. From the doc page:

"This FOH method differs from standard causal FOH ...."

The c2d implementation of foh looks like:

Hd = (1/T)*(z-1)^2 / z * Z(Linv(H(s)/s^2))

>> Hd=tf([1 -1],[1 0],T)*tf([1 -1],1,T)*c2d(zpk(H)*zpk(tf(1,[1 0 0])),T,'impulse')/T^2

Hd =
 
  0.04388 z (z-1)^2 (z-0.9048) (z+0.8752)
  ---------------------------------------
      z (z-1)^2 (z-0.8187) (z-0.7408)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

>> minreal(Hd)

ans =
 
  0.04388 (z-0.9048) (z+0.8752)
  -----------------------------
      (z-0.8187) (z-0.7408)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

>> c2d(zpk(H),T,'foh')

ans =
 
  0.04388 (z-0.9048) (z+0.8752)
  -----------------------------
      (z-0.8187) (z-0.7408)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

Keep in mind that the methods I'm showing using CST functions to implement the zoh, foh, and impulse conversions are most likely not how the CST actually implements c2d. I'm just trying to demonstrate the basic math.

Finally, I really like what you said: "Of course we do not just take this for granted."

[Dominik Hiltbrunner]

ZOH: It looks like you got the right answer for this case using the approach with the Symbolic Toolbox. You can also do the computation directly using Control System Toolbox functions:

I tried that too but I didn't manage to implement it using CST commands. So I went for a workaround using symbolic analysis. Also,

Hd=zpk(tf([1 -1],[1 0],T)*c2d(H*tf(1,[1 0]),T,'impulse')/T)

is a very elegnat method.

Impulse Invariance: You should add a comment to your table that the formula you show is only valid for the form of H(s) in the table, i.e., H(s) does not have repeated poles. Repeated poles have to be handled differently.

My text book doesn't handle cases with repititive poles. I guess "handled differently" refers to how poles with multiplicity > 1 are handled using partial fraction decomposition, i.e. there we chose a different approach too. However, it might be a much better idea to write the impulse invariant method differently, namely:

Because that way we can use this result to compute the zoh and foh (like you did).

Pole-Zero Matching: I don't think your formula is quite right. [...] you also need to include an additonal N-M-1 zeros at z = -1

You're right. My formula is not correct. Adding the following lines in my code fixes the problem.

for n=1:(N-M-1)
    Num = Num * (z+1);
end

And correcting the formula:

The reason is that your original transfer function has N-M zeros at s = j*inf, so you're discrete time transfer function should have the same number of zeros at the highest frequency point on the unit circle, which is z = -1.

Can you explain that? Why has the transfer function in the s-domain zeros at infinity on the imaginary axes? Soundly like I'm missing some fundamental property of transfer functions where .

[Dominik Hiltbrunner]

Furthermore, here is the updated table:

[Paul]

Regarding zeros at infinity for the pole-zero matching case ....

The zeros of the continuous time transfer function are the values of s where the numerator equals zero. However, if N > M, then for large s H(s) = 1/s^(N-M) (asuming leading coefficient of denominator is one) and H(jw) tends to zero as w -> inf. It is in this sense we mean that H(s) has N-M zeros at j*inf. Of course, the same would be true if you evaluate H(s) in any direction to infinity, not just along the jw axis. The difference N-M is important in that it governs the asymptotic slope of the gain curve of H(jw) and the asymptotic phase of H(jw) as w -> inf. When discretizing, the entire jw axis maps to the upper half of the unit circle. So in that sense the point z = - 1 is the highest frquency (at least the highest that we care about). So if you want the frequency response of Hd(z) to match that of Hc(s) at high frequency, one might argue that we want the same number of zeros at infinity of Hd(z) as there are in H(s), i.e., N-M zeros at z = -1. However, that would result in Hd(z) having relative degree zero (numerator and denominator having the same order), which means the input of Hd(z) has a direct feed through to the output, i.e., y(k) depends on u(k). But such a relationship may not be realizable due to computational delays in the processor, so the next best thing is to add N-M-1 zeros at z = -1 so that the Hd(z) has relative degree of 1, i.e., y(k) only depends on values of the input up to u(k-1), and lose a little bit of fidelity in the discretization. For example:

>> Hc = zpk(-1,-[2:5],1)

Hc =
 
           (s+1)
  -----------------------
  (s+2) (s+3) (s+4) (s+5)
 
Continuous-time zero/pole/gain model.

>> Hd=c2d(H,T,'matched')

Hd =
 
         0.00013342 (z-0.9048) (z+1)^2
  -------------------------------------------
  (z-0.8187) (z-0.7408) (z-0.6703) (z-0.6065)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

>> Hd1 = Hd*zpk(-1,[],1,T)/2  % divide by 2 to account for the extra gain at z = 1 due to the additional zero

Hd1 =
 
         6.6711e-05 (z-0.9048) (z+1)^3
  -------------------------------------------
  (z-0.8187) (z-0.7408) (z-0.6703) (z-0.6065)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model

>> bode(Hc,Hd,Hd1),grid
>> legend('Hc','Hd','Hd1')

You should see on the Bode plot that Hd1 does better at matching the phase of Hc(s) as w approaches the Nyquist frequency.

[Rik]

Hi Paul, sometimes the spam filter is a bit too strict, sometimes it misses spam. These things happen. I have now removed it from the spam filter.

[Dominik Hiltbrunner]

Hello Paul

Thank you for your detailed explainations and your effort. I hope I can return the favor one day.

I think all open questions are now answered.

[Tamas Kis]

I've made a quick function that performs the forward and backward transformations since they aren't supported by MATLAB:

https://www.mathworks.com/matlabcentral/fileexchange/90267-euler-c2d-transformation-c2d_euler