GeezerVison is showing github.com — fitfrd function · Issue #550 · python-control/python-control
DescriptionHi,
I am a student in learning speech signal processing. The control package helps me a lot. When I use python to reproduce a function, a problem bothers me for a long time.
I want to ask if there is a function can fit frequency response data with state-space model in python-control, which is similar to the fitfrd function in matlab. It would be too kind if you can provide some direction.
Sincerely. Activitymurrayrm commented on Feb 22, 2021 That functionality is not currently present, but it looks like it is implemented in the SLICOT SB10YD function, so it would be relatively easy to add. If you have some time and want to work on this, that would be great! Otherwise, we can mark this as "help wanted" and see if anyone has time to work on it. JinmingChe commented on Feb 22, 2021 Thank you very much for your reply! I will follow your advice to try to implenment the fitfrd function. If I succeed, I am happy to share my results. Before that, friends who can work on it are welcome. JinmingChe commented on Feb 23, 2021 Sorry for bother again. When I implenment the fitfrd function, using the control.freqresp function, it said i should use an interpolating FRD if I want additional points. But it doesn't show how to change the norm frd to interpolating FRD. I checked the document and there is no relevant explanation. murrayrm commented on Feb 23, 2021 Looking at the code, it looks like if you create the FRD object using smooth=True, the it defines an interpolating FRD. JinmingChe commented on Feb 24, 2021 I import the slycot package. In the dir of slycot, I can't find the SB10YD function. It's so kind to tell me how can find the SB10YD function. murrayrm commented on Feb 24, 2021 You'll likely need to write a wrapper around SB10YD in the slycot repository and generate a pull request there. JinmingChe commented on Mar 8, 2021 Hi, I find the sb10yd.f in slycot repository, and try to write a wrapper around SB10YD.
I would like to ask if I can wrap the sb10yd.f if I define the sb10yd function to the analysis.py and _wrapper.py files. In addition, under what circumstances should the parameter need to add the hidden?
If the question is too simple, please forgive me for being a newer. murrayrm commented on Mar 8, 2021 I think that is all that is required. But I've not written slycot wrappers before, so you may need to open an issue in that repository to get an answer. JinmingChe commented on Mar 8, 2021 Many thanks for your reply. bnavigator commented on Mar 8, 2021 I replied to python-control/Slycot#153 just before reading this thread. Please feel free to make an attempt and we can guide you along. The wrappers for SB functions should go into synthesis.pyf and synthesis.py. JinmingChe commented on Mar 15, 2021 Hellow, when I use the control.freqresp function it gives me three output (mag, phase and omega). But the freqresp funciton in matlab only one output in plural form. If there has function that I can transfer the frequency response output form ？ murrayrm commented on Mar 15, 2021 The following will work:
H = ct.rss(2, 1, 1)
omega = np.logspace(-2, 2, 100)
response = H(omega * 1j)
JinmingChe commented on Mar 15, 2021 My understanding is that you calculated the response from a state space object. When I run the code, the error message is displayed: TypeError:'StateSpace' object is not callable. murrayrm commented on Mar 15, 2021 What version are you running? This should work in the latest release (0.8.4). In that release and later, all LTI systems are callable and return the value of the transfer function. JinmingChe commented on Mar 15, 2021 My control version is 0.8.4. It also shows TypeError: 'StateSpace' object is not callable.
I'm not sure what you mean is that I can directly use frd sys to calculate the frequency response of the corresponding omega. Or I still need to use freqresp function to calculate mag, phase, omega, and then convert them to plural form. murrayrm commented on Mar 15, 2021 Oops. I was wrong: this is something that will be in version 0.9.0, which we will release in about a week.
I checked and the matlab module also returns mag, phase, freq instead of the complex value, which makes it inconsistent with MATLAB's usage. I'll open a separate issue on this. murrayrm commented on Mar 16, 2021 Correction to my prior note about fresp: this has been deprecated in v0.9.0 and so the right way to get this information is to call the LTI system with the complex frequency (or list of frequencies) that you want evaluated. See issue #573 for thoughts on adding fresp back into the MATLAB compatibility module (and making the output values consistent with MATLAB usage). JinmingChe commented on Mar 29, 2021 Hi, Matlab implements the .’ operation on the state-space model. How can I implement it in the control?
example (matlab):
sys = ss(Amat,Bmat,Cmat,Dmat,Ts);
sys = sys.'
JinmingChe commented on Mar 29, 2021 And I have another question. When I use the tf to create a transfer function system and then use the ss create a state space system. Compared with matlab, the control results are different.
example：
matlab —— ss(tf(n,d,-1))
result
A =
x1 x2 x3 x4 x5 x6
x1 -1.093 0.4225 0.8036 0.1752 -0.263 -0.2787
x2 2 0 0 0 0 0
x3 0 1 0 0 0 0
x4 0 0 1 0 0 0
x5 0 0 0 1 0 0
x6 0 0 0 0 0.5 0
B =
u1
x1 0.25
x2 0
x3 0
x4 0
x5 0
x6 0
C =
x1 x2 x3 x4 x5 x6
y1 0.1349 0.0666 0.03996 -0.0869 -0.1135 0.00342
D =
u1
y1 0.0851
control —— control.matlab.ss(control.matlab.tf(n.flatten(), d.flatten(), -1))
result
A = [[-1.09275252e+00 8.45135768e-01 1.60749491e+00 -3.50284601e-01
5.26406269e-01 2.78501895e-01]
[ 1.00000000e+00 -1.88597374e-17 4.11457136e-16 9.42161929e-16
-5.25031488e-17 2.81947516e-16]
[ 0.00000000e+00 1.00000000e+00 -6.29307309e-17 4.71371808e-16
7.54550775e-16 8.18501034e-17]
[ 0.00000000e+00 0.00000000e+00 -1.00000000e+00 8.37026245e-16
-5.02862888e-16 9.14800774e-16]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00
1.43285427e-16 5.61733088e-16]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
1.00000000e+00 -5.55111512e-16]]
B = [[1.]
[0.]
[0.]
[0.]
[0.]
[0.]]
C = [[ 0.03369084 0.03331414 0.01999275 0.04343376 0.056765 -0.00086861]]
D = [[0.08511411]]which n,d is the same as matlab. bnavigator commented on Mar 29, 2021
Compared with matlab, the control results are different.
State space realizations of transfer functions are not unique.
#431 (comment) bnavigator commented on Mar 29, 2021
Hi, Matlab implements the .’ operation on the state-space model. How can I implement it in the control?
You could Implement StateSpace.transpose() and make StateSpace.T an alias for it.
You would also have to make a distinction (and document!) between what MATLAB does for sys' vs sys.'. AFAICT the resulting systems both have the inputs and outputs exchanged, but sys' reverses the phase, while sys.' does not. JinmingChe commented on Mar 29, 2021 Thank you very much for your advise, it solved my problem very well. Another problem trapped me. When I use the hinfsyn function, the internal command runs to sb10ad (out = sb10ad(n, m, np_, ncon, nmeas, gamma, P.A, P.B, P.C, P.D)), debugging cannot be performed, the console does not display content, and there is no error message.
using command control.hinfsyn(P,1,1).
P is calculated using the control.augw() command.
slycot.version'0.4.0.0'
control.version'0.9.0'
Pycharm in window system. bnavigator commented on Mar 29, 2021 When it comes to usage, independent of implementing a fitfrd function, please use the "Discussions" forum.
Another problem trapped me. When I use the hinfsyn function,
--> #597 roryyorke commented on Jul 14, 2024 roryyorke commented on Jul 20, 2024 here's a naive first go at fitfrd:
import control as ct
from slycot import sb10yd
def fitfrd(gfrd, n, *, tol=0, minphase=False):
discfl = 0
flag = int(minphase) #
rfrdat = gfrd.fresp[0,0].real
ifrdat = gfrd.fresp[0,0].imag
omega = gfrd.omega
lendat = len(gfrd.omega)
print(discfl, flag, lendat, rfrdat, ifrdat, omega, n, tol)
nout, a, b, c, d = sb10yd(discfl, flag, lendat, rfrdat, ifrdat, omega, n, tol)
if nout != n:
# really raise a warning; can't trigger this
raise ValueError(f"asked for order {n}, got {nout}")
return ct.ss(a, b, c, d)With this fourth-order system with a notch, it doesn't do too badly.
SB10YD (or our interface?) doesn't handle DC gain < 0. I also couldn't easily get non-minimum-phase results.
I'm also not sure how to trigger nout != n.
import numpy as np
import control as ct
import matplotlib.pyplot as plt
from importlib import reload
import fitfrd
reload(fitfrd)
from fitfrd import fitfrd
from control import StateSpace
from numpy import array
g = ct.tf([1], [1, 0.1, 1]) * ct.tf([1], [1, 0.1, 16]) * ct.tf([1,0.1,4],[1])
states = 4
pz = np.concatenate((abs(g.poles()),
abs(g.zeros())))
eps = 1e-7
if max(pz) == 0:
raise RuntimeError('seems implausible')
wmin = min(pz)/10
wmax = max(pz)*10
omega = np.geomspace(wmin, wmax, 1 + 20 * states)
gfrd = ct.frd(g, omega)
h = fitfrd(gfrd, states)
hmp = fitfrd(gfrd, states, minphase=True)
labels=['g','h','hmp']
print(f'{any(h.zeros().real > 0)=}')
plt.clf()
mags, phases = ct.bode([g, h, hmp], omega=omega).flat
for mag, phase, label in zip(mags, phases, labels):
mag.set_linewidth(3)
mag.set_alpha(0.5)
mag.set_label(label)
KybernetikJo commented on Jul 22, 2024
I'm also not sure how to trigger nout != n.
Returns
-------
n : int
The order of the obtained system. The value of n
could only be modified if n > 0 and flag = 1.So basically only for a minimal phase system n can be chosen.
SB10YD (or our interface?) doesn't handle DC gain < 0
You can also compare the behavior with
octave control:fitfrd
which also uses sb10yd. roryyorke commented on Aug 4, 2024
You can also compare the behavior with
Good idea, thanks. I installed Octave 9.2.0 and octave-control 4.0.1, and it has much the same behaviour, i.e., low-frequency gain assumed to be positive. (One could test 1/s vs -1/s, 1/s^2 vs -1/s^2, etc., but I think this is sufficient to demonstrate the issue.)
version = 9.2.0
==========
1/(s+2) LF response (freq, resp):
0.1000 + 0i 0.4988 - 0.0249i
fitfrd result for 1/(s+2):
Transfer function '' from input 'u1' to output ...
0.006587 s + 0.9988
y1: -------------------
s + 1.997
Continuous-time model.
==========
-1/(s+2) LF response (freq, resp):
0.1000 + 0i -0.4988 + 0.0249i
fitfrd result for -1/(s+2):
Transfer function '' from input 'u1' to output ...
0.006587 s + 0.9988
y1: -------------------
s + 1.997
Continuous-time model.
octave:3> pkg list control
Package Name | Version | Installation directory
--------------+---------+-----------------------
control *| 4.0.1 | /home/rory/.local/share/octave/api-v59/packages/control-4.0.1
whoops, edit to add the script:
% Demonstrate LF DC gain bug in fitfrd, probably due to SB10YD
display(version);
pkg load control
%%%% 1/(s+2)
fprintf('\n==========\n\n');
gtf = tf(1,[1,2]);
gfrd = frd(gtf);
[h, w] = frdata(gfrd);
fprintf('1/(s+2) LF response (freq, resp):\n');
display([w(1), h(1)]);
fprintf('fitfrd result for 1/(s+2):\n');
display(tf(fitfrd(gfrd, 1)));
%%%% -1/(s+2)
fprintf('\n==========\n\n');
gtf2 = tf(-1,[1,2]);
gfrd2 = frd(gtf2);
[h2, w2] = frdata(gfrd2);
printf('-1/(s+2) LF response (freq, resp):\n');
display([w2(1), h2(1)]);
fprintf('fitfrd result for -1/(s+2):\n');
display(tf(fitfrd(gfrd2, 1)));
MetadataMetadataDevelopmentNo branches or pull requestsIssue actions
|
