Confirmation on how to calculate phase differences at every frequency point between two wideband signals?

Question

I have a problem that I thought was going to be simple, but it has become surprisingly stubborn and I question my method...

I have used the method (described below), but I wanted to confirm that this was the right approach, and see if there was another method I could use.

What I basically have is two wideband signals, (on the order of 100 Mhz BW), in quadrature. So I have the I and Q of those two signals. What I am trying to find is the phase difference between those two signals as a function of frequency.

The method I am currently trying is the following:

1. I take my quadrature complex signals, and window each of them with a hamming window 3 times over.
2. I take the FFT of the windowed complex signals, at some arbitrary length, say 2048. Call the results $S_1(w)$ and $S_2(w)$.
3. I multiply $S_1(w)$ by the conjugate of $S_2(w)$.
4. I take the argument out of the complex exponential vector from the above step (inverse tangent of imag over real). This gives me the phase differences between the two signals as a function of frequency.

My questions on this method:

• Is this the right approach for this simple type of problem?
• When I tried this method for a simulation with 1, 2 or a couple tones, I can get good results. I find the index of where my tones exist, (look at the abs of the FFT), I take that index, and look at the $\Delta$phase value from step 4 above. This gives me the right phase.
  • However, when I look the $\Delta$phase values for frequency indicies where said those frequencies dont exist, I get wild answers. Why is this? I realize those frequencies dont exist but shouldn't the phases just be 0 - 0 = 0 in that case?
• Based on this (with a couple tones), can I readily extend this to a wideband signal and also expect good results?
• My last question regards the nature of the phase delay here - physically what is happening is that one of the signals (optical) is going through a device that delays different frequencies differently. Of course the delays are manifested as phase offsets between all frequencies... but at the same time there is a time delay associated with this as well. Even if we were to measure the $\Delta$phase between the two signals properly, wouldnt it still be impossible to measure the time delay between them at those frequencies? I can have a time delay of 10 wavelenghts, but a phase offset of 0.

Thoughts?

Thanks in advance!

EDIT FROM COMMENT FEEDBACKS:

I think I am clearer now as to what my question really is: For a case where you have well separated tones and high SNR, my method can work.

However, what to do in the wideband case? In this sense, suppose $s_1[n]$ and $s_2[n]$ both contain two tones that are very close to each other in the frequency space. But your FFT can only give you so much frequency resolution. So instead of applying my method to two peaks, you only have one peak to go on. Now what?

• What relation would the $\Delta$phase of that peak have to the actual $\Delta$phase of the 'true' peaks? Is it their average?
• One method I am thinking of is correlating the two signals in the time domain, and working backwards to get the phase offsets per frequency taking advantage of the fourier translation property, although part of me thinks it might not work since there is no new information, just looking at the problem from a different domain.
• Another thought I had for the wideband case would be to put the whole thing through many narrow band filters, and then apply my original method. Is there a name for something like this?

Thanks again

Answer 1

What you seem to want is the phase response of the device that your optical system is passing through. Based on that, your high-level approach makes intuitive sense overall; I'm not sure why you chose a Hamming window and applied it three times. In any event, though, if you have a linear discrete-time system with (complex-valued) frequency response $H(e^{j\omega})$, then:

$$ Y(e^{j\omega}) = H(e^{j\omega})X(e^{j\omega}) $$

where $X(e^{j\omega})$ and $Y(e^{j\omega})$ are the discrete-time Fourier transforms (DTFTs) of the input and output signals, respectively. If you express the above in polar format (magnitude and phase), then the phase response of the system is just:

$$ \angle H(e^{j\omega}) = \angle Y(e^{j\omega}) - \angle X(e^{j\omega}) $$

that is, the phase response of the system just expresses the amount of phase shift that the system applies to the input signal as a function of (normalized) frequency $\omega$. Your approach takes advantage of the fact that for two complex numbers $X(e^{j\omega})$ and $Y(e^{j\omega})$:

$$ \angle Y(e^{j\omega}) - \angle X(e^{j\omega}) = \angle \left(X^*(e^{j\omega})Y(e^{j\omega})\right) $$

So, given the DTFTs of the system's input and output signals, you can straightforwardly calculate the phase response that you seek. The devil is in the details of how one would obtain the DTFT representations of these signals. The nice part about the DTFT is that it is continuous in frequency; the bad part is that the DTFT sum is infinite in length. So, if you really want to evaluate the phase response at arbitrarily many frequencies, you better be patient, because you'll need to collect an infinite amount of data first.

The discrete Fourier transform (DFT) is often used as an approximation to the DTFT, which is valid under certain conditions. If the signals $x[n]$ and $y[n]$ are finite in length (which for a practical scenario like yours, they are), then the DFT can be thought of as evaluation of the DTFT on the finitely-spaced frequency grid $\frac{2\pi k}{N}, k = 0, 1, \ldots , N-1$ (where $N$ is the signal length). So, the larger you make $N$, the smaller your sample spacing across the frequency axis.

Most physical devices will have a frequency response that is smooth to some degree (taken with the caveat that I'm not knowledgeable on optical device physics), so you should be able to find a DFT length that gives a decent approximation to the DTFTs for the input and output signals of your system. You can then apply your approach to estimate the phase response of the device.

---

You noted in your question that you see incorrect results for the estimated phase response for DFT bins in the output signal that didn't contain any content in the input signal; this is to be expected. Your assumption that you would just get zero phase for those bins is not likely to hold in a practical scenario; noise present in the sampled signals, spectral leakage, or other interference can cause you to observe some small amount of energy at frequencies that you wouldn't expect it. These extraneous sources can introduce distortion between the input and output that is not due to the device's response, and for complex numbers with small magnitude, small amounts of additive noise can drastically affect the number's resulting phase.

Therefore, you can only consider the phase measurement valid at frequencies where there is significant spectral content in the input. Stated differently, if you want to make a precise phase measurement, you want high signal-to-noise ratio in the complex value that you're using to estimate the phase.

Answer 2

Noise pretty much has a random phase. So it's usually reasonable to just clamp the phase of any bin with a magnitude below some noise threshold to zero. Make sure this threshold is above the (unavoidable) numerical noise level.