# Obtaining the open loop gain - (estiamte the gain and phase frequency response)

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

1. We turn off the feedback loop. We apply a 'tone' or some type of controlled input signal.
2. We calculate the power spectral densities of the input and output signals
3. We calculate the cross-spectral density between the input and output signals.
4. We then calculate the ratio of the cross-spectral density to the power spectral density of the input signal
5. 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)

• According to this, csd has two return values. I am not familiar with python but what happens when you write pxy = csd(x,y) instead of f, pxy = csd(x,y) ?
– AJN
Aug 17, 2023 at 2:36
• If you don't get an answer here, you can try asking at the signal processing SE.
– AJN
Aug 17, 2023 at 2:37
• Yes, csd does return 2 values. One is the range of frequencies, the other is the cross-spectral power density at each frequency bin. The line, pxy = pxy[1,:], selects all of the cross-spectral densities. If I wanted, frequencies = pxy[0,:], but this will just be 0 to f_sample. Aug 17, 2023 at 12:25
• Based on the post here, dsp.stackexchange.com/questions/89065/…, I incorporate the change needed. Assuming that the way I am estimating the transfer function is correct, are the steps listed above sufficient for obtaining the open loop gain? Aug 17, 2023 at 14:36
• Welch's method implementations usually also allows one to input window size, window function, and amount of overlap between consecutive windows. Does the Python implementation you are using also supports this? It might also be interesting to look at the coherence between x and y. Aug 19, 2023 at 0:33