Finding phase of a received signal using spectrogram function

Question

I am trying to find the phase of my dataset between a range of frequencies. I can't post my code or data but generally speaking, the frequencies range from 10 to 100kHz and I need the phase between 20 and 30kHz.

The data is read into Python as real numbers in the form of 64-bit floats. I was thinking to first pass the data through a spectrogram and set the mode to angle or phase. My first question is, I see from the scipy documentation the difference between the angle and phase modes is the wrapping. I'm not sure which to use so what do I need to consider when choosing one?

Second, how should I choose the parameters of the spectrogram to get the correct phase i.e. how should I choose the overlap number, fft length and window?

Lastly, once I have the complex angle returned by the spectrogram function, what operations do I need to do to recover the phase in degrees? Or, is there a better way to do this that does not use the spectrogram function?

Thank you

Answer 1

1. If complex-valued STFT is of interest, preferably (and in this case definitely) avoid scipy, librosa, MATLAB; the only suitable STFT implementation I know of is ssqueezepy's (disclaimer, am author).
2. Any hop_size > 1 will alias STFT, and unlike for a spectrogram, it's always deleterious for phase extraction, unless aliased phase is acceptable.
3. No one method can handle all signals, since time-frequency decomposition is non-unique, hence so is phase decomposition. If your signal is expected to be composed of sufficiently separated AM-FMs, then time-frequency methods (STFT, CWT) will work well, and an upgrade can be synchrosqueezing.

I'll go with assumption in 3. The task is analogously "intrinsic mode" extraction, i.e. finding $x_1$ and $x_2$ in $x(t) = x_1(t) + x_2(t)$. This is accomplished by "carving" time-frequency; what this means, and tools to achieve it, are described here. ssqueezepy doesn't have a built-in carving function for non-SSQ_CWT, so the one-integral inverse must be understood and tested on full-signal inverse if you are to write code for anything other than SSQ_CWT. For STFT it should suffice, however, to just correctly select the 2D coordinates.

Once inverses are obtained, it becomes standard phase recovery via the analytic signal - so in the referenced post, exclude the .real step, and apply angle and handle unwrapping (e.g. np.unwrap, but search around the site for more info to be sure).

Demo

Adding below code to ridge_chirp.py produces:

where hilbert would not be needed if issq_cwt didn't use .real. Note, "recovered" is from "original" plus WGN with variance 2.

# NOTE: also use `noise_var = 2`
from scipy.signal import hilbert

po, prec = [np.unwrap(np.angle(hilbert(g))) for g in (xo, xrec)]
pkw = dict(show=1, w=.6, h=.8, xlabel="time [sec]")
plot(ts, xo)
plot(ts, xrec, title="Signal: original vs recovered", **pkw)
plot(ts, po)
plot(ts, prec, title="Phase: Original vs recovered",  **pkw)