2
$\begingroup$

This is related to this post: https://engineering.stackexchange.com/questions/56050/obtaining-the-open-loop-gain-estiamte-the-gain-and-phase-frequency-response

Using python, I would like to write a method which estimates the transfer function similar to the MATLAB tfestimate(...) found here: https://www.mathworks.com/help/signal/ref/tfestimate.html#bvi03i4-1

The goal is to estimate the transfer function based only on observing the input and output sampled data. For the below, x is the sampled input data going into my 'black box' system. y is the sampled output data coming out of my 'black box' system. My target waveform are 8k in length, however the below snippet should work for any size as long as x and y are the same length (bigger should be better and the input should be a chirp or tone in the operational range).

Python code Snippit:

# Calculate the power spectral densities of the input and output signals.
    pxx = welch(x)
    pyy = welch(y)
    # Calculate the cross-spectral density between the input and output signals.
    pxy = csd(x, y)
    
    # Convert the cross-spectral density so that it plays nice with numpy
    pxy = np.asarray(pxy)
    
    #Convert input and output power spectral destines to np arrays so we can divide them laer
    pxx = np.asarray(pxx)
    pyy = np.asarray(pyy)
    f = pxx[0,:]
    pxx = pxx[1,:]
    pyy = pyy[1,:]
    pxy = pxy[1,:]
    
    # Calculate the ratio of the cross-spectral density to the power spectral density
    # of the input signal.
    h = np.divide(pxy, pxx)

Does this look correct?

Other helpful links that don't work for me but capture the idea:

https://stackoverflow.com/questions/28462144/python-version-of-matlab-signal-toolboxs-tfestimate

Transfer function estimation from frequency response

$\endgroup$
1
  • $\begingroup$ Did I answer your question? $\endgroup$
    – Jdip
    Aug 23, 2023 at 7:00

1 Answer 1

1
$\begingroup$

It looks almost correct. You might end up with a reversed phase.
The correct expression for $\tilde{H}(\omega)$, which estimates the transfer function $H(\omega)$ is: $$\tilde{H}(\omega) = \frac{P_{yx}(\omega)}{P_{x}(\omega)}$$

In code, you should have:

# Calculate the cross-spectral density between the input and output signals.
        pyx = csd(y, x)

See also this answer and this one for a lot more detail.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.