fitfrd function

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.

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.

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.

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.

You'll likely need to write a wrapper around SB10YD in the slycot repository and generate a pull request there.

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

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

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

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

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 fitfrd function, 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

octave control:fitfrd

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

Navigation Menu

Repository navigation

fitfrd function #550

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

fitfrd function #550

Description

Activity

murrayrm commented on Feb 22, 2021

JinmingChe commented on Feb 22, 2021

JinmingChe commented on Feb 23, 2021

murrayrm commented on Feb 23, 2021

JinmingChe commented on Feb 24, 2021

murrayrm commented on Feb 24, 2021

JinmingChe commented on Mar 8, 2021

murrayrm commented on Mar 8, 2021

JinmingChe commented on Mar 8, 2021

bnavigator commented on Mar 8, 2021

JinmingChe commented on Mar 15, 2021

murrayrm commented on Mar 15, 2021

JinmingChe commented on Mar 15, 2021

murrayrm commented on Mar 15, 2021

JinmingChe commented on Mar 15, 2021

murrayrm commented on Mar 15, 2021

murrayrm commented on Mar 16, 2021

JinmingChe commented on Mar 29, 2021

JinmingChe commented on Mar 29, 2021

bnavigator commented on Mar 29, 2021

bnavigator commented on Mar 29, 2021

JinmingChe commented on Mar 29, 2021

bnavigator commented on Mar 29, 2021

roryyorke commented on Jul 14, 2024

roryyorke commented on Jul 20, 2024

KybernetikJo commented on Jul 22, 2024

roryyorke commented on Aug 4, 2024

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions