This question is attempting to consolidate material discussed in two previous topics. These are:
Understanding the H1 and H2 estimators @Pontus S
I am estimating Complex Frequency-Domain Transfer Functions.
I will call the frequency-domain system X H = Y.
I am comparing estimators:
H0 = Y/X (i.e. direct frequency domain division OUTPUT/INPUT)
H1 = (YX*) / (XX*) (* is conjugate, H1 is supposedly best for output noise)
H2 = (YY*) / (XY*) (best for input noise)
Note that I only have a single recording of each experiment, so I cannot do any averaging of the spectra in H1, and H2.
For synthetic and real data, the three estimates appear identical. (Not surprising as they are mathematically identical, as pointed out by @Pontus in Q 71811.)
- Is the reason my H0, H1, H2 estimates are the same because I am not able to do any averaging in the spectra. Hence I am not reducing the influence of noise.
- In this situation is there anything to be gained by using H1, H2. rather than just H0 (Direct division)?
Feedback to @V.V.T.
Thanks very much for the effort you have put into this. It is appreciated. However, your code is investigating a different problem to the one in my post. You are feeding three different signals to the three estimators. My question relates to real data, where the noise is already added. If you modify your code slightly you will see what I mean. Here is what you get when you apply the three estimators all to same signal, with noise on the output. The three estimates plot on top of each other. i.e. they are identical
CODE: import numpy as np import scipy from scipy import signal import matplotlib.pyplot as plt from scipy.signal import max_len_seq #from scipy.fft import fft, ifft # 'INPUT' sequence x = max_len_seq(8)*2.0-1.0 plt.figure(0) plt.plot(x) plt.show() lx=len(x) # Comb filter b, a = signal.butter(3, 0.05) #awgn noise z = np.random.normal(0, 0.1, len(x)) y = signal.lfilter(b, a, x) # signal with noise on output############ yn = y + z plt.plot(yn) plt.show() #H1 estimator auto/cross spectral densities pxx pxy _, pxy = signal.csd(x,yn) _, pxx = signal.csd(x,x) H1 = np.divide(pxy,pxx) plt.figure(1) plt.title("Noise on Output") plt.plot(abs(H1), label="H1") #H2 Estimator f, pyy = signal.csd(yn,yn) f1, pyx = signal.csd(yn,x) H2 = np.divide(pyy,pyx) plt.plot(abs(H2), label="H2") #H0 estimator (Spectral division) flat=np.ones(lx,float) #just to enable consistent use of csd #to get the spectrum of x, and yn _, px = signal.csd(x,flat) _, pyn= signal.csd(yn,flat) H0 = np.divide(pyn,px) plt.plot(abs(H0), label="H0") plt.legend() plt.show() # signal with noise on input####################3 xn = x + z y2 = signal.lfilter(b, a, xn) _, pxy = signal.csd(xn,y2) _, pxx = signal.csd(xn,xn) #H1 Estimator H1 = np.divide(pxy,pxx) plt.figure(2) plt.title("Noise on Input") plt.plot(abs(H1), label="H1") #H2 estimator, auto/cross spectral densities pyy pyx f, pyy = signal.csd(y2,y2) f1, pyx = signal.csd(y2,xn) H2 = np.divide(pyy,pyx) plt.plot(abs(H2), label="H2") _, px = signal.csd(xn,flat) _, pyn= signal.csd(y2,flat) H0 = np.divide(pyn,px) plt.plot(abs(H0), label="H0") plt.legend() plt.show() exit()
I am not very familiar with the signal.csd routine, but it is possible that with the short sequence used here (255) there is no averaging/windowing being done. This is my situation (no spectral averaging). Note that when the signal is made longer, using the csd code some differences emerge, and this may be because spectral averaging is being done in csd. This is what what suggested by @Pontus in his final comment (#74836), and the purpose of my question was to see if he or anyone else could confirm this with real-data experience.
Again thanks very much for your time and effort. I have upticked your answer to reflect your input.