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
