How to properly smooth the phase of a spectrum (or any unit-complex function)

Question

I want to smooth the phase of a measured (transfer) spectrum without destroying unit-complexity of the phase factor. Suppose $$f:\mathbb{R}\to \mathbb{C}\qquad , \qquad f(\omega)=r(\omega)\cdot {\rm e}^{i\phi(\omega)}$$ is the spectrum, then it is well-known how to smooth the power spectrum (i.e. $f^*f=r^2$): just multiply the autocorrelation function (i.e. inverse Fourier transform of $f^*f$) with a window function $w(t)$ with limited width and use the convolution theorem to show that this is the same as filtering $f^*f$ with the Fourier transform of the window. For example, if $w(t)$ is a Gaussian, $w(\omega)$ will also be Gaussian (with reciprocal width with respect to the width of $w(t)$).

But what about smoothing the phase? If one simply smoothes the phase factor (which has unit modulus) in the same way, the result usually will not have unit modulus anymore. Moreover, phases that change rapidly but linearly will cause vast cancellation. But fast phases in a transfer function correspond to big translations/lags of wave groups, which I want to be represented faithfully in the transfer function. If I just cut them off (which is what happens if I filter the phase factor), bigger translations/lags will just disappear.

For example take $$f(\omega)=e^{i\omega t_0}$$ as a pure phase function, then averaging this over a frequency interval from $\omega=\omega_0-\pi/t_0\dots \omega_0+\pi/t_0$, the result will be zero (in complex) independent of $\omega_0$ because it always integrates over a whole unit circle: $$\bar{f}(\omega_0)=\frac{1}{2\pi/t_0}\int_{\omega_0-\pi/t_0}^{\omega_0+\pi/t_0} e^{i\omega t_0}d\omega = 0$$ By contrast, the average phase angle will be linearly increasing, i.e.: $$\bar{\phi}(\omega_0)=\frac{1}{2\pi/t_0}\int_{\omega_0-\pi/t_0}^{\omega_0+\pi/t_0} \omega t_0 d\omega = \omega_0 t_0 = \phi(\omega_0)$$

So what I actually want to smooth is the phase angle $\phi(\omega)$ itself! But joining phases continuously is very challenging for measurement data, so the ATAN2 function is not an option. What I have so far come up with is the following. Suppose we had a pure phase function $$f(\omega)={\rm e}^{i\phi(\omega)}$$ Then $$f^\prime(\omega) f^*(\omega)=i\phi^\prime(\omega)$$ So if $\phi(\omega)$ is a linearly increasing phase angle, $\phi^\prime(\omega)$ is constant, independent of how fast the angle changes. Hence, $\phi^\prime(\omega)$ is perfectly destined to be smoothed by convolution, in my opinion. I don't know yet how to properly reconstruct a (smoothed) phase factor from that, other than directly integrating it.

Does this make sense? I want make sure it does before I calculate any further. I haven't seen anything like that in the literature. Maybe this method is already known. Of course, what I would most enjoy in the end is a simple smoothing formula that does not depend on whether $f$ has unit modulus or not.

Edit: Upon further thinking, I have recognized that the computation of the phase gradient (aka group delay) via $f^\prime(\omega) f^*(\omega)=i\phi^\prime(\omega)$ is exact only for a continuous spectrum, while it can get very inaccurate when applied to a discrete spectrum. The derivative for the discrete spectrum becomes a finite difference and if the phase changes very rapidly, the angles between adjacent spectral points become very blunt, making the finite differences inaccurate. So there is no way around computing the phase explicitely and unwrapping it, like Dan Boschen has proposed in his solution.

Answer 1

Update: After chatting with @oliver this turned out to be more challenging than I anticipated.

My approach for comparison against improved alternatives is as follows:

Unwrap the phase to create samples of unwrapped phase versus frequency

Filter this vector with a moving average filter (directly in the frequency domain) or any other interpolation /smoothing appraoch. I used a zero-phase filter (filfilt()) specifically to eliminate having to compensate for what would be a frequency shift.

Example:

fdata= randn(1000,1) + j*randn(1000,1);   # original frequency data
y = unwrap(angle(fdata));
y2 = filtfilt(ones(10,1),10,y));          # zero phase 2 pass moving average
fdata2 = abs(fdata)exp(j*y2);             # smoothed data

Here is the result zoomed in showing the smoothing of the phase.

I would be interested in another answer that outlines a simple mathematical approach on the original complex data samples without converting to unrwapped phase vs frequency and back.

My (modified) first inclination was to divide the time domain samples in to consecutive blocks and sum the DFT of each block zero padded out to the original length : If the DFT can be evenly divided by M (such that N/M is an integer) then the sum of the M-length DFT's would result in exactly every Mth sample of the single DFT of the original sequence. Zero padding simply interpolates the samples in between. So this works as a frequency domain interpolation but cannot preserve the phase (except for the phase of every Mth sample) unless the group delay (which is the negative derivative of phase with respect to frequency) is constrained over M samples such that the deviation is constrained to $\pm \pi$ in that interval. (which is the frequency domain equivalent of memory over multiple samples in the time domain). if samples in that interval deviate more than $\pm \pi$ they cannot be uniquely distinguished from samples that deviated less than than. This is an interesting counter-part to white noise in the time domain, where white means the spectral density in the frequency domain is constant, where if the spectral density was constrained to low frequencies we could easily interpolate phase from sample to sample in the time domain. So here in this case if we have "white-noise" in the frequency domain samples-- this means the power density over time is flat.... but then it stands to reason that if we constrained the time domain to the low times (which means removing the delay due to causality as well which manifests as a linear phase slope in frequency) we should be able to resolve the interpolated phase values accurately.

Here is a plot comparing the unwrapped phase over all samples to that with every other sample from a randomly generated phase signal in frequency (each frequency bin is selected from a complex white Gaussian noise process - so uniform phase distribution). The x-axis is frequency.

In comparison here is a result for data that is windowed and time with all delay removed (again the x-axis is frequency):

Here is a plot of the magnitude of the time domain data for the second plot (the magnitude for the first plot would be evenly distributed, random) demonstrating the condition under which phase samples can be interpolated in frequency. This is a non-causal zero-delay waveform with limited time response. (A delay in time is a linear-phase slope in frequency and thus would contribute to the condition of modulo $2\pi$ transitions in phase between frequency samples which can then no longer be uniquely resolved):

Answer 2

You can take a wavelet transform of the given data and obtain a a phase angle time series and then average out the phase angle. The important thing to keep in mind is that this averaging is not the same as averaging a magnitude spectrum. To see an example have a look at the figure below:

So if you were to simply average the angle values ($2\pi$ and $-0.1$), then you would get the angle shown in brown vector which is incorrect. So while averaging we have to make sure that we average the vectors of those angles and thus obtain the yellow vector shown in the figure. Keeping this in mind, we obtain the phase angle time series from the wavelet decomposition. The morlet wavelet is suitable for this. Once we have the angle time series we average out the angles (moving average filter for example) using the concepts illustrated above.