How to Average Complex Responses (and Justification)? I am developing software that calculates the response of a system by comparing the FFT of input and output signals.  The input and output signals are divided into windows and, for each window, the signals are median-subtracted and multiplied by a Hann function.  The instrument response for that window is then the ratio of the FFTs of the processed data.

I believe the above is standard procedure, although I may be describing it poorly.  My problem comes in how to combine the responses from the multiple windows.

As far as I can see, the correct approach is to average the complex values, across all windows.  The amplitude and phase response are then the amplitude and phase of the average, complex value at each frequency:

av_response = sum_windows(response) / n
av_amplitude = sqrt(real(av_response)**2 + imag(av_response)**2)
av_phase = atan2(imag(av_response), real(av_response))


with implicit loops over frequency bins.

But I have been asked to change this to calculate amplitude and phase in each window first, and then average the amplitudes and phases across all windows:

amplitude = sqrt(real(response)**2 + imag(response)**2)
av_amplitude = sum_windows(amplitude) / n
phase = atan2(imag(response), real(response))
av_phase = sum_windows(phase) / n


I have argued that this is incorrect because averaging angles is "just wrong" - the average of 0 and 360 degrees is 180, for example, but the people I am working with responded by saying "OK, we will only display amplitude".

So my questions are:


Am I correct in thinking that the second approach is generally incorrect
for amplitudes too?
If so, are there any exceptions that may be relevant, and which may explain
why the people I am working with prefer the second method?  For example,
it looks like the two approaches will agree as the noise becomes small, so perhaps this
is an accepted approximation for low noise?
If the second approach is incorrect, are there any convincing, authoritative references that I can use 
to show this?
If the second approach is incorrect, are there any good, easy to understand examples that show this for amplitude
(as the average of 0 and 360 degrees does for phase)?
Alternatively, if I am incorrect, what would be a good book for me to educate
myself better?


I have tried to argue that the average of -1 1 1 -1 1 -1 -1 should be zero rather than 1, but that was unconvincing.  And while I think I could, with time, construct an argument based on max likelihood estimation given a particular noise model, it is not the kind of reasoning that the people I am working with will listen to.  So, if I am not wrong, I need either a powerful argument from authority or an "obvious" demonstration.

[I tried to add more tags, but can't find relevant ones and can't define new ones as a new user - sorry]
Transfer function estimation is usually implemented slightly differently than the method you describe.

Your method computes

$$\left\langle \frac{\mathcal{F}[y]}{\mathcal{F}[x]} \right\rangle$$

where $\langle$angle brackets$\rangle$ represent averages taken over data segments, and a windowing function is applied to each data segment before taking the Fourier transform ($\mathcal{F}$).

A more typical implementation will compute the cross spectral density of x and y divided by the power spectral density of x:

$$\frac{\langle \mathcal{F}[y] \cdot \mathcal{F}[x]^* \rangle}{\langle|\mathcal{F}[x]|^2\rangle} = \frac{\langle \mathcal{F}[y] \cdot \mathcal{F}[x]^* \rangle}{\langle\mathcal{F}[x]\cdot\mathcal{F}[x]^*\rangle}$$

Where $\cdot$ represents a pointwise product, and $*$ the complex conjugate.

I believe this is to reduce the effect of data segments where bins of $\mathcal{F}[x]$ are excessively small.

Incoherent estimation

Your employer has suggested that you estimate the transfer function using 

$$\frac{|\langle \mathcal{F}[y]|\rangle}{\langle |\mathcal{F}[x]|\rangle}$$

This will work, but has two big disadvantages:


You don't get any phase information.
If your measurements of the input $x$ and output $y$ have any additional noise, then the transfer function estimation will not be correct.


Your method and the method I described circumvent these problems by using coherent averaging.

References

The general idea of using overlapped, averaged segments to compute power spectral densities is known as Welch's method.  I believe the extension to using this to estimate transfer functions is also often known as Welch's method, although I'm not sure if it is mentioned in Welch's paper.  Looking up Welch's paper might be a valuable resource.  A useful monograph on the subject is Bendat and Piersol's book, Random Data: Analysis and Measurement Procedures.

Validation

To validate your software, I suggest applying several test cases, where you generate Gaussian white noise and feed it through a digital filter with a known transfer function.  Feed the inputs and outputs into your transfer function estimation routine and verify that the estimate converges to the known value of the transfer function.