How to properly smooth the phase of a spectrum (or any unit-complex function)
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):