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

  • $\begingroup$ 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. $\endgroup$ May 5, 2023 at 15:05

1 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).


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)
  • $\begingroup$ 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? $\endgroup$ May 4, 2023 at 17:28
  • $\begingroup$ Updated. No ifft is needed. $\endgroup$ May 5, 2023 at 14:47
  • $\begingroup$ 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). $\endgroup$ May 6, 2023 at 16:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.