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

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.

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):

• Although I have the feeling that this may probably be a part of the solution, I might point out that I haven't said anything about a discrete spectrum (which I will obtain during practical computation, though). If you could make this - somehow - work also for a continuous Fourier transform, it might be it. – oliver May 1 at 13:32
• I see-- I am confused as you mentioned measurement data so not really sure how it can be truly continuous ---so you are not using numerical techniques of any kind (which would all be discrete) and are looking for a formulaic solution? Also with processing data if you unwrap the phase you avoid any of the roll-over issues. – Dan Boschen May 1 at 13:36
• You are right, it is about measurements, but even if I start out with discrete data, padding zeros in the time domain in order to interpolate the spectrum gives me what? By stating the problem for the continuum, I don't have to care about apparent solutions that are just equivalent to resampling. I actually want to reduce detail in the spectrum, which I already do by smoothing the PSD, but I want it also for the phase. But maybe I just don't understand the consequences of my requirements well enough... – oliver May 1 at 13:44
• I will keep thinking about your suggestion, though... – oliver May 1 at 13:45
• I see you want to smooth the data you already have--- then this implies you are sampling higher than you need in which case I'll add more detail to that. – Dan Boschen May 1 at 13:46

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.

• I appreciate your unusual view to the problem. I have tried to tackle wavelets several times, but always ended up at "too complicated", "too unpredictable", "too many degrees of freedom", while I know exactly what to expect from Fourier transform. I think even if I understood more about wavelets, I would not know how to implement what you have suggested. Maybe you can elaborate some more... – oliver May 1 at 13:38
• Hnmm, the comments would fall short to explain all of that, let me give you a link to a similar video, he explains it in the context of obtaining consistency over trials, but you can apply the same concept to filter the phase angles, I have not found wavelet to be unpredictable, remember, FFT is good for stationary signals mostly and it's popularity is because there exist a fast algorithm available, the wavelet too has fast algorithms available and it has features that one simply cannot see in an FFT analysis, youtu.be/R1Pro555H6s – Dsp guy sam May 1 at 13:47
• I have reviewed the video (and others in this series) and it seems to me that what you/the video are saying is nothing but a restatement of my problem. Say we have a phase function $f(\omega)=e^{i\omega t_0}$, and we want to average this phase factor over a frequency interval from $\omega=\omega_0\dots \omega_0+2\pi/t_0$, the result will always be zero (in complex). Contrarily, the average phase angle will be linearly increasing, as is obvious from the definition of $f$. Hence, just averaging the 'phase vectors' does not give the right result nonetheless. – oliver May 1 at 18:07
• You are suppose to average accross time, it's the phase angle "time" series" – Dsp guy sam May 1 at 18:49
• I'd rather understand what I am supposed to do. – oliver May 1 at 20:03