Repository navigation
fitfrd function #550
Description
Activity
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.
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.
Reacted by Ben GreinerSorry 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.
Looking at the code, it looks like if you create the FRD object using
smooth=True, the it defines an interpolating FRD.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.
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.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.
Many thanks for your reply.
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.
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 ？
The following will work:
H = ct.rss(2, 1, 1) omega = np.logspace(-2, 2, 100) response = H(omega * 1j)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.
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.
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.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
matlabmodule 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.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 addingfrespback into the MATLAB compatibility module (and making the output values consistent with MATLAB usage).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.'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.0851control ——
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.
Compared with matlab, the control results are different.
State space realizations of transfer functions are not unique.
#431 (comment)Hi, Matlab implements the
.’operation on the state-space model. How can I implement it in the control?You could Implement
StateSpace.transpose()and makeStateSpace.Tan alias for it.You would also have to make a distinction (and document!) between what MATLAB does for
sys'vssys.'. AFAICT the resulting systems both have the inputs and outputs exchanged, butsys'reverses the phase, whilesys.'does not.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.When it comes to usage, independent of implementing a
fitfrdfunction, please use the "Discussions" forum.Another problem trapped me. When I use the hinfsyn function,
--> #597
Slycot has a SB10YD interface thanks to python-control/Slycot#203, released as part of Slycot v0.6.0.
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)
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
which also uses sb10yd.
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.1whoops, 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)));
Hi,
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.