# freqresp gives wrong output for purely real inputs

5 views (last 30 days)
Will Robertson on 17 May 2020
Edited: Paul on 4 Jun 2020
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)]
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 on 22 May 2020
Could you share the link for where you saw the comment "don't use freqresp -- evalfr gives proper results here"?
Will Robertson on 22 May 2020
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.

Devineni Aslesha on 26 May 2020
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.

Zhao Wang on 28 May 2020
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 on 28 May 2020
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 on 3 Jun 2020
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.

Paul on 24 May 2020
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 on 25 May 2020
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 on 25 May 2020
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 on 25 May 2020
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...