Recovering phase from spectrogram

Question

In this answer, an iterative algorithm is given for recovering a signal from its spectrogram, assuming overlapping windows:

$$ x_{n+1} = \operatorname{istft}(S \cdot \exp(i \cdot \operatorname{angle}(\operatorname{stft}(x_n)))), $$

where $S$ is the spectrogram and $\operatorname{(i)stft}$ is the (inverse) short-time Fourier transform. I tried this out using NumPy, SciPy, and librosa:

import librosa
import numpy as np
import scipy

def mse(phase_pred, phase_true):
    '''
    Calculate the mean square error between the true phase and the
    predicted (reconstructed) phase.
    '''
    return np.mean(np.angle(phase_pred/phase_true)**2)

# Load an audio file and calculate STFT.
x, sample_rate = librosa.load('audio.wav', sr=44100)
D = librosa.stft(x)
mag, actual_phase = librosa.magphase(D)

# Try to reconstruct the phase using the iterative algorithm above.
phase = np.exp(1.j * np.random.uniform(0., 2*np.pi, size=actual_phase.shape))
x_ = librosa.istft(mag * phase)
print('iter {} mse {}'.format(-1, mse(phase, actual_phase)))
for i in range(100+1):
    _, phase = librosa.magphase(librosa.stft(x_))
    x_ = librosa.istft(mag * phase)
    print('iter {} mse {}'.format(i, mse(phase, actual_phase)))
    if i % 10 == 0:
        scipy.io.wavfile.write('recons{:05d}.wav'.format(i), 44100, x_)

As far as I can tell, it isn't converging:

However, the audio after the 100th iteration definitely sounds closer to the original than the one with randomized phases.

Why isn't this algorithm converging to the correct phase? Am I misunderstanding its purpose?

Answer 1

It should be clear from the definition of the STFT that the result is invariant under global phase changes. With the restriction to real input signals that still provides for a sign ambiguity, or a phase difference of $\pi$. That means, even if your reconstruction is perfect, the result may differ by that phase from the original.

Another reason for the failure to converge may be the existence of nearly disconnected components in your signal. Say you have two sound components that are temporally separated so that there is no STFT frame that overlaps with both. Each component may be reconstructed with a different sign, giving you four different reconstructions with identical STFT. The same works for components that are separated well enough in frequency direction. Each component doubles the ambiguity of the reconstructed signal and will increase the probability of a deviation from the original signal.

The final reason I can see lies in the nature of the iterative algorithm. Fixed point iterations, of which kind this algorithm is, can fail to converge to a proper fixed point, keep oscillating or even diverge. Unless the maps involved are well behaved in a certain sense, predicting convergence is hard. It seems that your iteration doesn't settle, and it's unclear if adding iterations will help. So this may apply here, together with the other two problems.