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_overlap_percent - these will influence frequency resolution vs. PSD estimator variance trade-off.