This post is related to the following stackexchange link but in this new question, I am focusing on describing the controller and going over my process of obtaining the open loop vs closed loop gain: controls - how to measure phase margin for low gain system
Description of Plant: The plant is a black box which looks similar to a band pass filter. However, this system varies over time and is difficult to directly model. The goal is to apply and control a well regulated electric field.
Description of Controller: We are working with a simple PI controller which attempts to regulate an electric field's strength. This is a digital system and we have the ability to see the input into the controller and the resulting output. We can not continuously stream the input and output but we can look at 8,192 (8k) consecutive samples for both the input and output. This digital interface is the only way we have to measure what is going into or out of the system. The sample rate is very high (56 MHz) but in my previous post I normalized it (0 to 0.5 fs) to keep the plot legends readable.
My Crazy Science Project to obtain the Gain Response (estimate): The goal is to determine the gain/phase response by observing the actual input/output signals and estimate the transfer function. We are doing something similar to MATLAB tfestimate(x,y).
- We turn off the feedback loop. We apply a 'tone' or some type of controlled input signal.
- We calculate the power spectral densities of the input and output signals
- We calculate the cross-spectral density between the input and output signals.
- We then calculate the ratio of the cross-spectral density to the power spectral density of the input signal
- This gives us a transfer function which can be used to produce the magnitude and phase plots.
Given the above, does this sound reasonable for estimating the open loop gain/phase response?
Thank you for all of your help.
Here is the Python code. Note, x is the 8k input waveform and y is the 8k output waveform. # 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(y, x) #NOTE, the input names are swapped, this is to prevent phase reversal.
# 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)
csd
has two return values. I am not familiar with python but what happens when you writepxy = csd(x,y)
instead off, pxy = csd(x,y)
? $\endgroup$