averaging power spectrum from multiple signal of different length

Question

I have multiple physiological signals (EMG - muscular electrical activity) that are all of different length (one may even be 2/3 times longer than the smallest one). These signals are from 2 types of muscular activity, all sampled at 1000Hz. So, say I have:

• 10 EMG signals from activity A, with an average lenght of 1.5 sec, ranging from 0.5 to 2.5 sec.
• 10 EMG signals from activity B, with an average lenght of 10 sec, ranging from 8 to 12 sec.

What I want is to compare the average power spectrum of activity A and B to show if there are any differences. So, my plan is to compute all 10 power spectra of signals from activity A and average them. Then, do the same for all 10 signals from activity B. Then plot the 2 averaged spectra for comparison.

But, as all signals are of different lengths, the power spectra are of different lengths as well. So how do I average them ? So far, I saw few solutions, but I am not convinced which one is the best or maybe there are other ones.

• For each activity to compare do: segment all signals in equal-length segments. For each signal, compute the power spectrum for all segments and average them. This gives an average power spectrum for each signal. Then, average over all signals. From here, it seems that this method is appropriate when each signal is stationary over all segments. Which is not the case for my set of signals as the activities are performed in a dynamic environnement.

• For each activity to compare do: either zero-pad all signals to the longest size or crop all signals to the smallest size. This is to make all signals have equal length (I would prefer the first option as there is information I don't want to lose). Then, compute the power spectrum for all signals and average them. The problem is that zero-padding adds fluctuations in the power spectrum.

• For each activity to compare do: keep the size of all signals unchanged. Then compute their power spectrum and apply a moving average on them. Then apply interpolation to get each spectrum to equal lengths, and average them over all signals.

EDIT: Actually, the point to compare the spectrums is to show that both activities are different enough so that they could be differentiated in a real-time machine-learning strategy, using a sliding window of, say, 150 millisecond, were each window is an instance to classify. Ultimately, their differentiability have to be tested with an machine-learning algorithm directly, but I need to show first that the strategy is legit (i.e. that there is potential differentiability). Basically, I am doing an amplitude comparison between activities and, among other schemes, it includes spectrum comparison.

Here is an example of 2 signals from both activities and their corresponding spectrum. Also, for each signal, the spectrum is computed with scipy as follow:

xf = fft.fftfreq(len(signal), sampleRate)
yf = np.abs(fft.fft(signal, len(sig)))
PSD = yf * yf
plt.plot(xf, PSD)
plt.show()

Answer 1

To add to @Hilmar's answer, and address your comment on his answer:

Why is this enough to state that the signal is "reasonable stationary" ? doesn't it depend on the size of the frame? indeed, there are relatively quiet portion in the B activity signals. Too small frame would see drastic change, doesn't it ?

Welch's method uses frequency averaging to estimate the PSD, so you need to be averaging frames that have somewhat consistent frequency content. That's the definition of stationarity. One way to deal with quasi-stationarity is to use median averaging instead of mean averaging. That way if a frame has mostly silence, the effect on averaging will be minimized. Scipy's Welch has an option to do that.

Also, the signals from activity A are relatively small in duration. is this still ok to do Welch's method with signals that short ?

It depends on the frequency resolution you're happy with. 1s of data at 1000Hz is 1000 samples. Say you compute the PSD using frames of length 128, the resolution will be ~8Hz. If that's ok with you, then you should be fine. In general though, the longer the frames, and the higher the number of frames you average together, the better the estimate.

For your purposes, I'd say you should be ok with Welch. If you actually had to get a robust estimate at lower resolution, you'd need to move on to different PSD estimation methods.

Answer 2

Provided that the signals are reasonable stationary (i.e. they are either have the "activity" or not and don't change state in the middle of the file), the standard Algorithm to estimate the PSD is Welch's method. See for example https://ccrma.stanford.edu/~jos/sasp/Welch_s_Method.html

You basically chop the signal into overlapping frames, apply a window, perform an FFT and then average all frame spectra in energy. This allows you to estimate the PSD for each signal with the same frequency grid and resolution, regardless of length. The longer the signal, the more accurate the estimate becomes. You also get to pick the frequency resolution that's meaningful and relevant for your application and not just what happens to be the sample rate divided by the signal length.

Both Matlab and Python have pre-fab functions for that: pwelch()