Zero-order hold practical differentiator offsets magnitude

This is more about theory than MATLAB, but this seems like a good place to ask. If I have a practical differentiator, meaning

where

and discretize it with zero-order hold I notice the magnitude plot starts to offset as (see figure below.) If I repeat this with a first-order hold or bilinear transformation I get what I would have expected, the magnitudes track the pure differentiator until near the cutoff frequency and then flattens.

I thought this was interesting and wondered if someone could shed some light as to why this offset is happening for zero-order hold? A big thank you in advance.

Note: The figure shows cutoff frequencies as Hz in the legend but rad/s in the x-axis. I just wanted to make this quick without much formatting.

[Paul]

Hi Greg,

Are you asking to see the math that yields the shown result? Or are you asking for a conceptual explanation that is consistent with the math?

[Greg]

That's a good question/point you've made. I'm sure I could crank out the math and get some insight as to what tethers a FOH/bilinear transform of this sysem to the low frequency magnitudes, and very well may do that. Really any insight you're willing to provide I would be happy to hear. I'm a Mech. E. that's just getting started with these topics and found it interesting that for ZOH the magnitude is shifted based on higher cutoff frequencies.

Answer 1

[Paul]

Hi Greg,

I don't know if this will be the complete explanation your looking for, but here goes ....

Define the continuous-time transfer function

syms s z
syms tau Ts positive
syms n integer
H(s) = s/(tau*s + 1)

H(s) = 

The zero order hold approximation is:

Hzoh(z) = simplify((1-inv(z))*ztrans(ilaplace(H(s)/s,s,n*Ts)))

Hzoh(z) = 

Its frequency response is

syms omega real  % rad/sec
Hfzoh = Hzoh(exp(1j*omega*Ts))

Hfzoh = 

At low frequency, we can substitute a Taylor series for exp(1j*omega*Ts)

Hfzoh = simplify(subs(Hfzoh,exp(1j*omega*Ts),taylor(exp(1j*omega*Ts),omega,0,'Order',2)))

Hfzoh = 

At low frequency as omega -> 0, Hfzoh can be further approximated by

Hfzoh = Ts*omega*1j/tau/(1 -exp(-Ts/tau))

Hfzoh = 

which is a pure differentiator scaled by a gain

zohgain = Ts/tau/(1 - exp(-Ts/tau))

zohgain = 

The Tustin approximation is

Htustin(z) = simplify(H(2/Ts*(z-1)/(z+1)))

Htustin(z) = 

Its frequency response

Hftustin = Htustin(exp(1j*omega*Ts))

Hftustin = 

Approximate with the Taylor series

Hftustin = simplify(subs(Hftustin,exp(1j*omega*Ts),taylor(exp(1j*omega*Ts),omega,0,'Order',2)))

Hftustin = 

As omega -> 0, we have

Hftustin = 2*omega*1j/2

Hftustin = 

which is the pure diffentiator.

Verify the ZOH result numerically

wnyquist = 1000*pi;    % rad/sec
Fs = wnyquist*2/2/pi;  % Hz
Ts = 1/Fs;             % sec

fc = 200;
tau = 1/(2*pi*fc);
Gs = tf([1 0],[tau 1]);
Gzoh = c2d(Gs,Ts,'zoh');
zohgain = Ts/tau/(1 - exp(-Ts/tau));
Gtustin = c2d(Gs,Ts,'tustin');

figure
h = bodeplot(tf('s'),Gs,Gtustin,Gzoh/zohgain);  % divide Gzoh by the zohgain derived above
h.PhaseVisible = 'off';
xline(1/tau,'g')

I noticed that the parameters in the problem are such that Ts > tau. Typically, we'd want the sampling period, Ts, to b 5 - 10x smaller than tau. If we use the low end of that rule-of-thumb, the zohgain would be

Ts = tau/5;
zohgain = Ts/tau/(1 - exp(-Ts/tau))

zohgain = 1.1033

and the ZOH comparison would much better at low frequency

Gzoh = c2d(Gs,Ts,'zoh');
Gtustin = c2d(Gs,Ts,'tustin');

figure
h = bodeplot(tf('s'),Gs,Gtustin,Gzoh);  % don't divide Gzoh by the zohgain
h.PhaseVisible = 'off';

So that's the mathematical explanation, i.e., the ZOH approximation introduces a gain into Gz, and that gain gets closer to unity as Ts becomes smaller relative to tau.

From a physical perspective, the best I have for this particular problem is as follows.

The Tustin approximation is just a straight substitution for s in terms of z, where z is defined by a low order approximation to exp(-s*Ts). At low frequencies where that approximation is valid, we'd expect Gztustin to recover the same frequency response as Gs.

I'm less sure of the situation with the ZOH approximation, so take what follows with a grain of salt.

The ZOH approximation is based on taking differences between step responses of the plant induced by input samples that are held constant due to the D/A converter. The step response of the plant is

syms tau s t Ts
assume([tau Ts],'positive')
H(s) = s/(tau*s + 1);
hstep(t) = ilaplace(H(s)/s,s,t)

hstep(t) = 

At t = 0, the first sample of input to the plant is u_0 (which is held constant from 0 <= t < Ts) and the the output of the plant due to the first sample is

syms u_0
sympref('HeavisideAtOrigin',1);
y_0(t) = u_0*(hstep(t)*heaviside(t) - hstep(t-Ts)*heaviside(t-Ts))

y_0(t) = 

The output of the plant due to the second sample is

syms u_1
y_1(t) = u_1*(hstep(t-Ts)*heaviside(t-Ts) - hstep(t - 2*Ts)*heaviside(t - 2*Ts))

y_1(t) = 

The output of the plant at t = 0 is:

y_0(0)

ans = 

The output of the plant at t = Ts+ (i.e., after the input at t = Ts is applied) is the sum of the output due to the first sample and the second sample

y_0(Ts) + y_1(Ts)

ans = 

Ideally we'd want that second output sample to be just u_1/tau so that the difference between the successive output samples is proportional to the difference between the input samples, just like differentiation. That ideal situation will only be approached as Ts becomes small relative to tau and that exponential term approaches unity. In other words, the sample period, Ts, has to be small enough reatlive to tau so that response of the plant due to u0 has sufficiently decayed by the time the second smaple u_1 is applied.

I didn't think too much about the FOH case, but I suspect that it's taking advantage of the fact that the FOH approximation is based on the ramp response of the plant.

[Greg]

Amazing, Paul! This was a very thorough and thoughtful response which completely answers my question. Seeing the gain terms really spells it out, and I had a similar inclination about ZOH responding to steps vs FOH responding to ramps. Thank you so much!

[Paul]

You're quite welcome.

I'm still not happy with the second half of the analysis that looked at the time domain response (the stuff with u_0 and u_1) and how it relates to the steady state output to a sinusoidal input. I'll think about that some more and post back if I come up with something better.