Real Data Complex Transfer Function Using H0, H1, H2 Estimators

UPDATE

Python's signal.csd routine, as it is used in my code, don't work miracles, at least in the sense you would expect ("but it is possible that with the short sequence used here (255) there is no averaging/windowing being done" -- as if it does averaging/windowing for >255 sequences). From the code you can see that, when presented only the input and output data, the program has no means to separate the output into pure signal plus noise and can only do its best to recover the transfer function of the LTI system from its input and output data as if the output signal plus noise is a newly defined output signal and there is no noise in the system at all -- to solve the ill-posed problem of deconvolution.

We can re-formulate the problem: multiple realizations of the LTI system's output with noise at the output are given for an identical input data. Providing that there is no noise at input, what is the best estimate for the system transfer function? The code below solves the problem with 64 realizations given.

import numpy as np 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)[0]*2.0-1.0 _, pxx = signal.csd(x,x) # linear filter b, a = signal.butter(3, 0.05) # noiseless output y = signal.lfilter(b, a, x) z = np.random.normal(0, 0.1, len(x)) # add noise to output yn = y + z _, pxy = signal.csd(x,yn) H1 = np.divide(pxy,pxx) plt.figure(1) plt.title("H1: a noise realization added to output") plt.plot(H1)

# average output noise over ensemble yn = y + z/64 for iter in range(1,64): yn = yn + np.random.normal(0, 0.1, len(x))/64 # cross psd with averaged output over ensemble _, pxy = signal.csd(x,yn) #H1 estimator H1 = np.divide(pxy,pxx) plt.figure(2) plt.title("H1: yn averaged over sixty four realizations of noise at output") plt.plot(H1)

#spectral division X = fft(x) Y = fft(y) H0 = np.divide(Y, X) plt.figure(3) plt.plot(H0[0:len(H0)//2]) #draw plots plt.show()

And it is exactly what ZR Han posited in the answer to your question.

But pay attention to the conditional "as it is used in my code": you can employ the independence of noise component samples to avail and, to some degree, replace averaging over a signal ensemble by averaging over these i.i.d. variables of a unique signal realization. This is the ergodic theory in action: under certain conditions, averaging over ensemble is equivalent to averaging over time.

The scipy.signal.csd function has the optional 'window' parameter. By default this parameter is 'hann', the window length of 256, but I doubt that default parameters result in the noise handling in the way required for estimators, that is, averaging over independent realizations. However, there is a window value, 'dpss', which can be used to apply said 'ergodicity' to your avail. I would not bore you with a lecture on this approach in my rendering, you can read elsewhere about the technique called a "multitaper" method.

To my knowledge, there is no readily available scipy's implementation of multitaper. However, you can implement it for yourself.

Start with a simpler task of implementing the psd calculation with multitaper. For the toy problem you can take an unsolved task from the question How to accurately compute the Winger-Ville Distribution of an exponential. The OP was upset that he cannot replicate numerically the Wigner-Ville distribution. The principal reason was incorrect generation of the Wiener process realization, but in the end he required the precise coincidence including the distribution tails, and it is impossible with spectral division exactly for the same reason that you receive "the same result" for estimators: the averaging over ensemble is required. You can use the distribution tail suppression task to acquire skill in multitaper application. You can start with a multitaper implementation for PSD calculation -- for example, taken from github, as e.g. https://github.com/xaratustrah/multitaper or else -- and later to extend it or write your own implementation for cross channel spectral distribution calculation. Or maybe you would accommodate scipy's csd, using window='dpss', into your solution from the beginning. If you engage in this task, make known to dsp.SE users when you post your solution to github.