evalfr for discrete systems gives strange results

Hi! I've got some strange behaviour from the function evalfr. Basically I use evalfr to extract the frequency response of a system in the form of a complex number (which has a magnitude and phase).

When I do this for a continuous system everything works fine. If I then discretize the system I get results I cannot explain. The obtained response for the discrete system appears to vary with the sampling time. (I check/compare the results to the bode function). Can someone take a look/explain?

%%continuous system (works as expected)
clear
close all

sys2 = tf([4700 4393 3.245e08],[1 7.574 1.202e5 0 0]);
bode(sys2), grid on

feval = 50                              %evaluation point in Hz
res2 = evalfr(sys2,feval*2*pi*i);       %extract frequency response as a complex number
mag2 = mag2db(abs(res2))                %convert magnitude to db
pha2 = angle(res2)*180/pi               %obtain angle and convert to degrees
[mag2c,pha2c] = bode(sys2,feval*2*pi);  %use the bode command to check
mag2c = mag2db(mag2c)
pha2c
%MAG2 AND MAG2C SHOULD MATCH (same for phase)

%%Discretized system (doesnt work as expected)
close all

T = 1/1000; %sample time
sys3 = c2d(sys2,T); %discretize
bode(sys3), grid on

res3 = evalfr(sys3,feval*2*pi*i);
mag3 = mag2db(abs(res3))
pha3 = angle(res3)*180/pi
[mag3c,pha3c] = bode(sys3,feval*2*pi);
mag3c = mag2db(mag3c)
pha3c
%MAG3 AND MAG3C SHOULD MATCH (same for phase)

Answer 1

[Chad MacDonald]

The evalfr command allows you to evaluate the frequency response of a system at any point in the complex plane. However, bode evaluates points only on the:

• Imaginary axis (s = jw) for continuous-time systems.
• Unit circle (z = e^jwT) for discrete-time systems.

Therefore, the reason you are getting inconsistent results between bode and evalfr is that you are using the s-domain complex variable with your discrete-time system rather than the z-domain variable.

The following code shows how to get the same response at a given frequency using bode and evalfr:

%%Create discrete-time system.
T = 1/1000;
sys = tf([4700 4393 3.245e08],[1 7.574 1.202e5 0 0],T);

%%Frequency = 50 Hz
w = 50*2*pi;

%%bode result
[magB,phaseB] = bode(sys,w)

%%evalfr result for z = e^jwT
resultE = evalfr(sys,exp(j*w*T));
magE = abs(resultE)
phaseE = angle(resultE)*180/pi

This information is not clearly conveyed in the documentation. I've updated the enhancement request in our tracking system, and we will update the documentation in a future release.