Skip to content

Use slycot's tb05ad for faster and more accurate frequency response #173

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
murrayrm merged 2 commits into from
Jan 6, 2018
Merged

Use slycot's tb05ad for faster and more accurate frequency response #173

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
138c097
Use slycot's tb05ad for faster and more accurate frequency response e…
rabraker Jan 2, 2018
62a29d9
Change omega_s to omega in statesp.freqresp()
rabraker Jan 4, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Change omega_s to omega in statesp.freqresp()
  • Loading branch information
@rabraker
rabraker committed Jan 4, 2018
commit 62a29d953dbb71d78c3c8dbab8e2b8855406e287
22 changes: 11 additions & 11 deletions control/statesp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -386,7 +386,7 @@ def horner(self, s):
return array(resp)

# Method for generating the frequency response of the system
def freqresp(self, omega_s):
def freqresp(self, omega):
"""Evaluate the system's transfer func. at a list of freqs, omega.

mag, phase, omega = self.freqresp(omega)
Expand All @@ -402,7 +402,7 @@ def freqresp(self, omega_s):

Inputs:
------
omega_s: A list of frequencies in radians/sec at which the system
omega: A list of frequencies in radians/sec at which the system
should be evaluated. The list can be either a python list
or a numpy array and will be sorted before evaluation.

Expand All @@ -414,30 +414,30 @@ def freqresp(self, omega_s):
phase: The wrapped phase in radians of the system frequency
response.

omega_s: The list of sorted frequencies at which the response
omega: The list of sorted frequencies at which the response
was evaluated.

"""

# In case omega is passed in as a list, rather than a proper array.
omega_s = np.asarray(omega_s)
omega = np.asarray(omega)

numFreqs = len(omega_s)
numFreqs = len(omega)
Gfrf = np.empty((self.outputs, self.inputs, numFreqs),
dtype=np.complex128)

@murrayrm murrayrm Jan 4, 2018

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?).

@rabraker rabraker Jan 4, 2018

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


# Sort frequency and calculate complex frequencies on either imaginary
# axis (continuous time) or unit circle (discrete time).
omega_s.sort()
omega.sort()
if isdtime(self, strict=True):
dt = timebase(self)
cmplx_freqs = exp(1.j * omega_s * dt)
if ((omega_s * dt).any() > pi):
cmplx_freqs = exp(1.j * omega * dt)
if ((omega * dt).any() > pi):
warn_message = ("evalfr: frequency evaluation"
" above Nyquist frequency")
warnings.warn(warn_message)
else:
cmplx_freqs = omega_s * 1.j
cmplx_freqs = omega * 1.j

# Do the frequency response evaluation. Use TB05AD from Slycot
# if it's available, otherwise use the built-in horners function.
Expand Down Expand Up @@ -476,8 +476,8 @@ def freqresp(self, omega_s):
for kk, cmplx_freqs_kk in enumerate(cmplx_freqs):
Gfrf[:, :, kk] = self.horner(cmplx_freqs_kk)

# mag phase omega_s
return np.abs(Gfrf), np.angle(Gfrf), omega_s
# mag phase omega
return np.abs(Gfrf), np.angle(Gfrf), omega

# Compute poles and zeros
def pole(self):
Expand Down