Transfer Function Estimate

Question

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

Answer 1

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.