0
$\begingroup$

I have been using scipy's spectrogram function (STFT) to compute a PSD spectrogram for a non-stationary signal. I provide X length window (with nperseg = X/10) each time to the spectrogram function and eventually concatenate the result. Now, I would like to use Welch for generating this spectrogram. I'm a bit confused on how to approach it.

Is there some working code that does it already?

$\endgroup$
2
$\begingroup$

Reading the documentation for scipy.signal.spectrogram I noticed that it does not do any kind of periodogram averaging. It simply splits up the signal into (possibly overlapping) segments, computes the magnitude of the DFT and plots it in each column of the STFT matrix.

For your application I would recommend looping through non-overlapping segments of the original (nonstationary) signal and using scipy.signal.welch to compute the Welch PSD estimate for each segment. For displaying the final result, you can concatenate these to form a matrix, where individual columns are the Welch PSD arrays for each segment. Here is some (mostly complete but untested) code to give you an idea:

# original non-stationary signal is in the array x
N = len(x)

# parameters to play with:
segment_length = 1024
segment_overlap_percent = 0
welch_window = 'hanning'
welch_nperseg = 256
welch_overlap_percent = 50

segment_jump_by = int(segment_length*(1.0-segment_overlap_percent/100.0))
welch_noverlap = int(welch_nperseg*welch_overlap_percent/100.0)

welch_specgram = np.zeros(#allocate space here)
for i in range(0,N-segment_length,segment_jump_by):
    current_segment = x[i:i+segment_length]
    f, welch_specgram_column = welch(current_segment, fs, nperseg=welch_nperseg, noverlap=welch_noverlap)
    # insert welch_specgram_column into welch_specgram at appropriate col number

#image plot welch_specgram matrix

There are some free parameters that you'll have to hand tune based on what trade-offs you are willing to make:

  • Length of each non-overlapping segment segment_length: longer segments will give you better frequency resolution in subsequent Welch periodogram estimation step but you don't want to make it too long otherwise you won't capture the non-stationarity of your signal.
  • welch_window, welch_nperseg, welch_overlap_percent - these will influence frequency resolution vs. PSD estimator variance trade-off.
$\endgroup$
  • $\begingroup$ Thanks. I will give it a try. I was just wondering, do you think that simple concatenation will work? Question is if the spectrogram will be displayed properly when each segment is computed and then concatenated individually or maybe I should have some overlap between different segments. $\endgroup$ – mike Oct 26 '17 at 8:53
  • $\begingroup$ If you use overlapping segments you'll introduce some correlation between their PSD estimates. If you are ok with that you can set segment_overlap_percent to something other than 0. See my edited code. $\endgroup$ – Atul Ingle Oct 26 '17 at 15:13
  • $\begingroup$ I was also wondering about another thing. The processing is performed in real time so every 10 seconds, I process a window by calling scipy.signal.spectrogram, but it returns "t" which is less than 10 seconds long. I would like to have this spectrogram in sync with some other data. Does it makes sense? Is there some way to deal with it? $\endgroup$ – mike Oct 31 '17 at 16:16
  • $\begingroup$ Not sure. My guess is that the time axis of spectrogram() may change slightly because of edge effects that depend on the spectrogram window length and overlap. $\endgroup$ – Atul Ingle Oct 31 '17 at 16:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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