Use slycot's tb05ad for faster and more accurate frequency response #173
Conversation
…valuation of state-space systems. If slycot is unavailable, use the built in horner function (instead of converting to a transfer function, as was done before).
control/statesp.py
Outdated
| @@ -386,22 +386,98 @@ def horner(self, s): | |||
| return array(resp) | |||
|
|
|||
| # Method for generating the frequency response of the system | |||
| def freqresp(self, omega): | |||
| """Evaluate the system's transfer func. at a list of ang. frequencies. | |||
| def freqresp(self, omega_s): | |||
There was a problem hiding this comment.
Is there a reason to use omega_s here instead of omega? In all other docstrings for freqresp (eg, xferfcn.py, lti.py) we just use omega. Also, shows up as omega below. I would leave as omega for consistency.
control/statesp.py
Outdated
|
|
||
| Inputs: | ||
| ------ | ||
| omega_s: A list of frequencies in radians/sec at which the system |
control/statesp.py
Outdated
|
|
||
| numFreqs = len(omega_s) | ||
| Gfrf = np.empty((self.outputs, self.inputs, numFreqs), | ||
| dtype=np.complex128) |
There was a problem hiding this comment.
Why complex128? I don't have a specific objection, but seems a bit arbitrary. OK to leave, unless there is some justification for doing something else (perhaps just `complex', which presumably defaults to machine precision?).
There was a problem hiding this comment.
I think I did this because the numpy page about dtypes doesn't list complex as a type.
Looking further down the page though, they say that dtype=complex will be interpreted as dtype=np.complex_ which is shorthand for np.complex128. To me, using np.complex128 is the most explicit as to what's happening.
This pull request addresses issue #116.
It uses the (relatively) newly wrapped TB05AD function from Slycot, when Slycot is available. If slycot is not available, I instead use the state space built-in method horner, rather than converting to a transfer function for the evaluation.