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

조회 수: 28(최근 30일)
Joana
Joana 2018년 11월 10일
편집: dpb 2018년 11월 29일
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
dpb 2018년 11월 11일
편집: dpb 2018년 11월 11일
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.
  댓글 수: 13
dpb
dpb 2018년 11월 29일
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.

댓글을 달려면 로그인하십시오.

추가 답변(0개)

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by