# Finding phase of a received signal using spectrogram function

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

• I also assume you confused "spectrogram" with "STFT", is that correct? spectrogram = abs(stft). If not, the problem is much more challenging. If yes, then "spectrogram" should be changed to "STFT" everywhere in the question. Commented May 5, 2023 at 15:05

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)

• Okay, let me know if I understand correctly. You are suggesting first to find the components of the signal using "carving" time-frequency which requires the one integral inverse but to use that I need the inverse FFT of the signal? From there I can then just apply the $angle$ and $unwrap$ functions since I will have complex data? Commented May 4, 2023 at 17:28
• Updated. No ifft is needed. Commented May 5, 2023 at 14:47
• I can show how to carve STFT in a separate Q&A, if you show an attempt of trying for yourself ("ordering code" is off-topic). Commented May 6, 2023 at 16:57