How to calculate PSD estimate of EMG data in 1Hz frequency bins for 100ms time windows.?

hi,

I have EEG signals of 1sec, and i want to calculate PSD in 1hz frequency bins for 100ms time windows. any advice will be appreciated :)

[dpb]

What was the sampling rate?

[Joana]

256

Answer 1

[dpb]

For those sampling parameters, one has that

Fs = 256;        % sampling frequency, samples/sec
T = 1;           % sample time, sec
dt = 1/Fs;       % sample interval, sec

From those we can derive that--

>> L=(T+dt)/dt         % length of overall time series
L =
 257
>> f=Fs*(0:(fix(L/2)))/L;  % frequency components of PSD (one-sided)
>> Fmax=f(end)             % maximum frequency of computed FFT, Hz
Fmax =
  127.5019
>> df=f(2)-f(1)            % frequency resolution, Hz
df =
    0.9961
>>

So, you have sufficient data to have approximately 1 Hz resolution in the full-sample PSD, but that encompasses the whole 1 second data acquisition period to get that much resolution.

If you limit to 100 msec, then you have only

>> L=round(0.1/dt)           % number of samples in approximately 100 msec real time
L =
    26
>>

Repeating the above with that for L and other sampling parameters as were, we find that

>> df=f(2)-f(1)       % frequency resolution of 100 msec time sample @ 256 Hz sample rate
df =
  9.8462
>>

this is, of course, about a tenth the resolution as it's one-tenth the previous time within the rounding of even/odd number of samples.

To get the 1 Hz frequency resolution over 100 msec intervals, you've got two choices --

1. Increase the sampling frequency by 10X
2. Interpolate by computing the N-based FFT instead

The second would be instead of the subsection of the signal as above

Y=fft(S(1:L));       % for first 100 msec, similar for subsequent

use

N=256;
Y=fft(S(1:L),N);     % compute N-length FFT on L samples of signal

I did an example using the example signal at doc fft without noise to illustrate what happens for the three cases; the figure is attached:

What is most apparent are the significant side lobes introduced by the interpolation process; windowing may help some, but it's an artifact of reality that you can't create information from nothing simply to be able to show a higher resolution on a figure.

While that is the most glaring artifact that shows up, the more significant result is that notice while there are 10X more points on the curve for the interpolated version, the shape of the peak itself is actually essentially the same as that of the short time sequence -- it illustrates you're not actually creating information and getting more resolution in the computation but for the most part simply interpolating between the points that are real data.

The one benefit is that you'll notice the peak for the 120 Hz component is closer to the actual frequency owing to there being frequency bins in the computed FFT that are closer to the true signal frequency. Note that in the original full-length that the power is split almost equally between two bins such that its magnitude must be estimated by the integral of the energy in the peak, not just with point value. In the example in the documentation, a signal length of 1500 was used instead of the 257 I used to match up with your sampling parameters to illustrate your case specifically.

To get the desired frequency resolution accurately you'll need to sample at a higher rate.

[dpb]

Added the original example results with Fs=1000 samp/sec, L=1500 samples so T=1500/1000 sec - 1dt = 1.499 sec. df=0.6667 Hz. The much better resolution owing to longer sample time is clearly obvious as is the precise magnitude of the estimate as numbers were chosen specifically so the bin spacing would exactly match the known frequencies and there is no noise.

[Joana]

Hi Thanks a lot for the detailed answer. I highly appreciate it.

If i have the epoched data (attached .mat file) , with fs=1200, in the order of: channel X Samples X Epochs, then how i can calculate the PSD for it, for each epoch. Because now i won't need the interpolation, right.?

[dpb]

Actually, if by "epoch" you're meaning the 100 msec windows, then the base frequency resolution is still 10 Hz because if you work through the algebra of the sampling relationships above

Fmax = 1/(2dt)  % Nyquist frequency, max in the signal

df = Fmax/(N/2) % frequency resolution

From latter, substituting in for Fmax,

df = 1/(N dt) = 1/T

The higher sampling rate will provide more samples so that the interpolation will be somewhat better, but the basic definitive resolution is limited by the length of the sampled signal, unfortunately.

There are some other techniques ("zoom" FFT) but the base DFT implemented in Matlab doesn't incorporate that.

As far as the mechanics of computing what you can with the base Matlab, if you reshape the data by L by column, then FFT will compute all the epochs in one call by column.

[Joana]

No by Epoch i meant 390 epochs. each epoch of 1second.

And isn't PSD estimate supposed to be by pwelch function or square of fft.?

Because in above comment you mentioned L point's fft, as a PSD calculation.?

[Joana]

and in here Y=fft(S(1:L),N);

S is the input signal, right.?

[dpb]

"390 epochs. each epoch of 1second."

Then where is the 100 msec coming from? I'm cornfoozed! :)

A 1 sec window will provide 1 Hz frequency resolution from baseband.

Welch or abs(fft) are two computational techniques to estimate PSD; sampling relationships are the same either way.

Yes, in the above I used S as the input signal, L is the actual signal length/number of samples and N is the FFT length; if choose to interpolate, N > L (in the olden days, one often used N=npow2(L) simply because compute power was so limited that there was a significant speed benefit in using the optimized power-of-2 algorithms instead of general. Now, with the available compute power that's essentially been relegated to history except for the most time-critical realtime situations or other very special circumstances).

[Joana]

10 windows each of 100ms=1sec

and then there are 390 epochs of 1sec

So i have to calculate PSD of 10x390 windows(of 100ms), in 1Hz frequncy bins.

[dpb]

Catenate the "windows" and reshape(S,[],390), then abs(FFT) (or PWELCH) will compute estimates by column.

[Joana]

Okay so here's the last query :)

i have fs=1200;

Time window= 250ms with 100ms overlap

so i will have L=300,

and if i don't go for the interpolation of increasing the sampling freq, i get the

Freq resolution= 4hz.

Right.?

Secondly of i use pwelch function like this : [pxx,f] = pwelch(tempData,window,overlap,nfft,sampligRate);

[pxx,f] = pwelch(Data,250,100,nfft,1200);

what will be the right choice for nfft.? 1024.?

[Joana]

I have written this code for calculation of PSD in 250ms time window, overlapped for 100ms.

Data is in this order: Number of channels X Number of samples X Number of iterations.

Now i have 60 channels, and 100 iterations, so what would be the order of pxx.? As with this code what i get is 513x1 (which 1024/2+1), shouldn't it be 60x513.?

sampligRate = 1200; % in 1 second
window = 250;
overlap = 100;
nfft = 1024;% 
psdBatch = [];
for epoch = 1:size(Data,3)
   for channel = 1:size(Data,1)
        % select samples to be analized
        tempData = Data(channel,:,epoch);
        % calculate PSD with welch method for specific channel and for specific
        %epoch
        [pxx,f] = pwelch(tempData,window,overlap,nfft,sampligRate);
        % output f shows the freq resolution
        psdBatch(channel,:,epoch) = pxx;
    end
end

[dpb]

Not clear to me what your data storage is, sorry...Attach a small file that illustrates with, say, three (3) channels and two (2) epochs. The magnitude of sizes isn't critical, it's understanding just what you've actually done in creating the data.

It would be more conducive to rapid calculations if you were to store channels column-wise with time the first dimension--that's the way all the ML functions are set up to work in a vectorized fashion.

[dpb]

Answer moved to comment -- dpb

<https://www.mathworks.com/matlabcentral/answers/uploaded_files/142152/Data.mat>

Hi,

Here's the data for 5 channels, and 100 epochs, wach of 1se (1200 samples). I need to calculate the PSD fetaures basically for all the channels and epochs.

please suggest me the better way to do it.

Because i am not sure if the code i have written is fine.

Also, i am not certain what should be nfft.? 1024.?

[dpb]

Well, again, we still aren't clear just what it is you really want.

What does ". I need to calculate the PSD fetaures basically for all the channels and epochs." mean in terms of the number of PSDs you would expect to compute? Taken literally, that would be 500, 5 channels x 100 "epochs", whatever those may be which is still not at all clear.(*)

Or, are we to be averaging over some number of those?

(*) Altho, if i do

D=permute(Data([2,1,3]));             % rearrange column by channel
plot([D(:,:,1);D(:,:,2);D(:,:,3)])    % plot "epochs" as extended times

it appears that perhaps each is a subsequent time fragment?

You can efficiently compute those FFTs from D as simply

Y=fft(D);

as FFT() is vectorized to treat a multidimensional array as individual columns down the first non-zero dimension. The result of Y is then 5x100 1200-element each complex double-sided FFTs just as we illustrated for a single vector earlier; you can convert to the one-sided PSD by the identical machinations on each column of 1200 elements and gather and average those as desired by plane or groups of planes ("epochs").

PWELCH() is NOT vectorized in the same way if you choose to use it instead; it can operate on 2D arrays by column so you'd have to iterate over the planes (epochs); again after permuting the storage as shown.

As for NFFT; that's your call; we talked quite extensively about the effects of interpolating or not.