freqresp gives wrong output for purely real inputs

I am stumped by the behaviour of evalfr vs freqresp in the control systems toolbox. They claim to give the same results, but when their arguments are purely real there is a (large!) discrepancy. I am inclined to believe the results from evalfr; for all sane transfer functions a real value of s should result in a real value of the transfer function.

With the following code:

G = tf([1],[1 1]); 
wc = 3+1e-10*1i; 
[evalfr(G,wc);freqresp(G,wc)]

Both give the answer 0.2500-0.0000i.

When the input is purely real, however:

G = tf([1],[1 1]);
wc = 3; 
[evalfr(G,wc);freqresp(G,wc)]

ans =
   0.2500 + 0.0000i
   0.1000 - 0.3000i

I'm using Matlab 2020a and I believe this behaviour has been around for a while; I came across an old code comment that said something like "don't use freqresp -- evalfr gives proper results here". Fine for me, but harder to advise when teaching students...

[Devineni Aslesha]

Could you share the link for where you saw the comment "don't use freqresp -- evalfr gives proper results here"?

[Will Robertson]

Sorry, I was unclear. That was my own code comment written to myself a few years back :)

FWIW the correct answer is definitely 0.25 here — it's literally as simple as 1/(s+1) for s=3.

Answer 1

[Devineni Aslesha]

Hi Will,

I have brought the issue of 'freqresp accepting frequency input with non-zero imaginary part' to the notice of our developers. They will investigate the matter further.

[Devineni Aslesha]

To compute the frequency response at purely real input frequency ‘wc’ using evalfr, we should pass the input frequency as ‘j*wc’ which is mentioned in the ‘Evaluate Frequency Response of Identified Model at Given Frequency’ example in evalfr documentation.

[evalfr(G,j*3) freqresp(G,3)] both gives the same results.

Also, when using freqresp, the frequency vector ‘w’ should be real to get the frequency response on the real frequency grid. If the frequency vector ‘w’ is complex-valued, its entries are interpreted as s or z values due to which [freqresp(G,3) freqresp(G,j*3) ] gives the same results.

[Will Robertson]

Yes, this behaviour has been established already in the earlier discussion on the page. There are two or three problems here:

1. freqresp is only documented to accept real values but it accepts complex ones anyway
2. freqresp is inconsistent in the results it returns since lim c->0 freqresp(r+c*1j) ≠ freqresp(r) — rather confusing
3. evalfr is described in the documentation as "a simplified version of FREQRESP" but this is misleading since the form of its inputs is different ( eval(G,3i)=freqresp(G,3) ) as it is in some sense generalised to complex inputs.

[Zhao Wang]

I would like to comment about the second point. If c is nonzero, no matter how small it is, the w is treated as a complex number. As a result, it is interpreted as s or z value in the transfer functions. The frequency response result would be different because the real value r would still be treated as the frequency value in freqresp(G,r).

You will get the same result using freqresp(G, r+0*1j) and freqresp(G,r) because r+0*1j and r are both treated as real values. For example, here is what I got:

>> freqresp(G,3)

ans = 0.1000 - 0.3000i

>> freqresp(G,3+0*j)

ans = 0.1000 - 0.3000i

[Will Robertson]

Thanks for the reply. I completely understand what is happening here (now). My main point in all of this now that my original confusion was cleared up, is that it just doesn't make sense to have a function that behaves this way. It is documented as accepting real inputs. It shouldn't accept a restricted set of complex inputs. If it needs to accept complex inputs for backwards compatibility, well, so be it :(

Your point about r+0i being real is completely correct from the Matlab perspective but is still confusing since my understanding is that the reals are a subset of the complex numbers. freqresp takes two "domains": reals, or "complex-but-not-reals". The latter domain is not a normal thing.

Maybe I should paint a picture here instead of using mathematical arguments. I had a problem I was solving from the literature that required evaluating a transfer function at desired locations for pole placement (therefore complex s). Maybe it was a copy paste error, maybe it was just a slip of the keyboard, but I ended up using freqresp to evaluate G at a few points. The frequencies were complex so with some quick sanity checks it was clear the code was working properly.

Then I tried placing a pole on the negative real axis. Results became instantly garbage. Switched to `evalfr`; code worked. Without stopping to check the docs (my bad) it was pretty obvious there was an inconsistency here and I raced to the discussion board to bring it up. Hopefully this discussion will be helpful to people who stumble across the same bump in the road...

As I said above, if backwards compatiblity were not an issue I would suggest only allowing real inputs into freqresp.

[Paul]

1. I agree with Will on point 1. The documenation is crystal clear. The input w must be real and the output is defined using the standard definition of freuqency response for either a continous or discrete system. I'm not sure why there's any reason for freqresp to accept an input with non-zero imaginary part. One might think that it should for backwards compatibility with the old version of freqresp (which is still usable, see my example below). But IMO that argument doesn't work because freqresp as currently implemented already has a different behavior than the old version for a real input. So that backwards compatibilty was lost anyway. At the very least, if freqresp is going to continue to accept complex inputs, it should generate a warning that that's unodocumented behavior, and it should be deprecated. And the same goes for the old version of freqresp.
2. Zhao, I can't comment on Will's point. But freqresp does return inconsistent results. Please check my example below that shows freqresp (the current version) returns different outputs for the same input. Surely that can't be acceptable.
3. As to Will's point 3, I agree that the documentation for evalfr is misleading. In fact, I think the name of the function is misleading. It should be evaltf. The very description says: "evaluates the transfer function ..." (emphasis added).
4. As has been well established we have for continuous time systems: evalfr(G,1j*w) == freqresp(G,w) with w real. I'd like to add that if G is a model for a discrete time system then we have: evalfr(G,exp(1j*w*G.Ts) == freqresp(G,w) if G.Ts is defined and w is real, and evalfr(G,exp(1j*w) == freqresp(G,w) if G.Ts is undefined.
5. The doc pages for evalfr and freqresp are pretty weak at explaining how these functions relate to each other, and the "More About" and "Algorithm" sections of doc freqresp are weak at explaining the discrete time case.

Answer 2

[Paul]

evalfr evalutes the transfer function at the value of the input argument, which can be an arbitrary complex number. For your example, evalfr(G,wc) = G(wc).

freqresp evaluates the frequency response at the value of the input argument, which must be real (doc freqresp). Because it's the frequency response, freqresp(G,wc) = G(1j*wc) (as long as wc is real).

So if you want to compute the frequency response of G using evalfr, you have to multiply the (real) input by 1j:

>> wc = 3; [evalfr(G,1j*wc); freqresp(G,wc)]

ans =

0.1000 - 0.3000i

0.1000 - 0.3000i

If you want to evaluate the transfer function of G at a complex number with non-zero real part, then your only choice is evalfr(G,s).

Check the doc pages for both functions.

Having said that, in 2019a the doc page for evalfr is missing the sections for the details on the input and output arguments. But it's all pretty clear from what is shown on that doc page.

[Will Robertson]

Hi Paul,

Thanks for your reply. The penny has now dropped — I knew that evalfr required the complex frequency input explicitly and at some point in my life I also knew freqresp was supposed to take a purely real input.

Where you say "If you want to evaluate the transfer function of G at a complex number with non-zero real part" that's not actually true:

w = 3+1i; [evalfr(tf(1,[1 1]),w); freqresp(tf(1,[1 1]),w)]
ans =

          0.235294117647059 -    0.0588235294117647i
          0.235294117647059 -    0.0588235294117647i

This example threw me off course and is arguably a bug in freqresp (that it allows a complex input at all) given its documentation states the input frequency should be real.

As currently implemented evalfr == freqresp for all inputs that are complex. For purely real inputs freqresp "flips" them to be complex. So to avoid confusion, like you say, don't use freqresp unless you are only looking for the response along the imaginary axis.

Regards,

Will

[Paul]

Will,

My statement that you quoted was intendended to be general in that I was thinking that 3, as in your example, is a complex number. I'm also suprised that freqresp (for a tf object like G, and probaby for its superclasses) accepts a frequency input with a non-zero imaginary part. And the situation is even worse than I thought. Consider:

>> G=tf(1,[1 1]);h1=freqresp(G,3);h2=freqresp(G,[3; 1+3*1i]);[h1 h2(1)]

ans =

0.1000 - 0.3000i 0.2500 + 0.0000i

freqresp computed a different output for the same input. Which is troubling. Though I can't say it's a bug because that second call to compute h2 is invalid anyway.

Then on top of this, the old-style use of freqresp is still usable, though I'm not sure if it's still supported functionality:

>> freqresp(1,[1 1],3)

ans =

0.2500

Apprarently freqresp used to work like evalfr works now.

[Will Robertson]

Sorry I misunderstood your statement about 3 being complex!

I guess it's a tricky to balance sensible behaviour with backwards compatibility. In a pragmatic world we probably can't change anything to avoid breaking old code...