computing delta power of EEG signal
显示 更早的评论
I wish to compute the “delta power” of EEG recorded over the course of a night’s sleep. From the raw EEG signal (in microvolts), I’m attempting to do the short-time Fourier transform (STFT) in small windows of the raw signal, then analyze the outputs in the range of 0.5-4Hz (delta).
I have lots of questions about 1) how STFTs should be ideally computed; 2) how Matlab computes STFTs; and 3) how to average results across frequency bands and time ranges of interest.
I apologize for the length of this post. It is clear that I should take a course in signal processing, but life is short, I’m in the middle of my career, and there is no time to learn things from 1st principles. Any help relieving my ignorance about all this will be greatly appreciated.
As background, I am running Matlab 7.5.0.342 with the Signal Processing Toolbox on a WinXP machine (64 bit).
1) My understanding is that I could take a spectrogram of my raw EEG signal as a way to implement STFTs. Does this seem reasonable? As opposed to snipping out small time windows of the raw signal, hanning windowing them, computing FFT, calculating power of output, moving on to next window, lather rinse repeat, etc…
2) For some reason, I have two spectrogram functions: “spectrogram” and “specgram”. The latter may be 3rd party. When I try to call the spectrogram, it gives me the error: Undefined function or method 'gencoswin' in ‘hamming’ at 14. Is one of specgram / spectrogram not included in the signal processing toolbox? Is this a licensing error? Is there any great difference between these two functions?
3) The few papers I’ve read on the process of extracting delta power breaks the signal into 4s windows (=nfft), computes the STFT, then averages across several of those windows (e.g., to get avg power in 20s windows). Some authors also average across frequency bin in order to get a single estimate of power in the range 0.5-4Hz. It is usually not specified how this averaging is done: since abs, square and log are non-linear, it matters when I average across frequency and time bins. Any suggestions at what step averaging should be done (e.g., should I average the log PSD or take the log PSD of the average?)
4) Related to 3), is the output b from specgram (or s in spectrogram) in squared units?
5) Right now, I’m using a window of the size of nfft. In this case, I assume there can be no overlap, is that true? Is there a problem not overlapping the signal?
6) Let's say I take a 20s snippet of my raw signal. Then I calculate a spectrogram with a window size (nfft) of 4s. I would expect there to be 5 time bins and Nyquist*4 (+1) frequency bins. I do get this result (e.g., f = Nyquist*4+1; t = 5), but when I plot the output (as surf), it goes from 0-16 sec, not 0-20 sec. That is, it is plotting 4 time windows from 0-4,4-8,8-12,12-16s. If I increase my time resolution (but worsen my spectral resolution), I get something closer to 20s windows. Is there an ideal trade-off? If I go with nfft=4s, then should I start the next 20s time window one 4s window behind, in order to make up for the lost window from 16-20s?
7) In older papers, the window is a boxcar the length of nfft. I think specgram takes a hanning window. Is there a way to specify boxcar for specgram?
8) Finally, another approach would be to do all this in the time domain: generate a, say, Chebyshev Type I bandpass filter between 0.5 and 4Hz; filter the raw EEG signal with these coeffs (using filtfilt to avoid phase shift errors); then simply compute 1og power on the filtered values of the form 10*log10(x2.^2 + 1) (the +1 is to remove fractions in the signal, which will produce negative values after logging). Should this produce the same thing as the spectrogram version above?
Thanks for the help
回答（8 个）
Kosai
2012-3-7
Hi Michael, I am going to answer the part of your question about the calculation of Delta-Power : at first you don't need to use the STFT to extract the Delta-Band-Information, but simply you have to use the Wavelet-Decompostion-Analysis(Look at the Wavelet-ToolBox), the most important point by wavelet analysis is to determine the Sampling-Freqeny(Fs) in your EEG-Data-Row, i will assume that your (Fs = 1000Hz as example) , so we need then an Wavelete-Decompostion-Analysis with (8 Levels) to achive the Band-Frequency = 4Hz(Delta-Band), simply type these commands :
S = "your EEG-Data-Row";
waveletFunction = 'db8' OR 'sym8' OR 'coif5' OR 'db4';
[C,L] = wavedec(S,8,waveletFunction);
%%Calculation The Coificients Vectors
cD1 = detcoef(C,L,1); %NOISY
cD2 = detcoef(C,L,2); %NOISY
cD3 = detcoef(C,L,3); %NOISY
cD4 = detcoef(C,L,4); %NOISY
cD5 = detcoef(C,L,5); %GAMA
cD6 = detcoef(C,L,6); %BETA
cD7 = detcoef(C,L,7); %ALPHA
cD8 = detcoef(C,L,8); %THETA
cA8 = appcoef(C,L,waveletFunction,8); %DELTA
%%%%Calculation the Details Vectors
D1 = wrcoef('d',C,L,waveletFunction,1); %NOISY
D2 = wrcoef('d',C,L,waveletFunction,2); %NOISY
D3 = wrcoef('d',C,L,waveletFunction,3); %NOISY
D4 = wrcoef('d',C,L,waveletFunction,4); %NOISY
D5 = wrcoef('d',C,L,waveletFunction,5); %GAMMA
D6 = wrcoef('d',C,L,waveletFunction,6); %BETA
D7 = wrcoef('d',C,L,waveletFunction,7); %ALPHA
D8 = wrcoef('d',C,L,waveletFunction,8); %THETA
A8 = wrcoef('a',C,L,waveletFunction,8); %DELTA
POWER_DELTA = (sum(A8.^2))/length(A8);
I hope, thats will help you ...
10 个评论
Kosai
2012-3-7
sorry but (Row = *Raw) :)
Michael
2012-3-8
Kosai
2012-3-8
hmm ok, in all cases i will not advice you to use the STFT Method to extract the Informations from the EEG-Data bcz you are going to lose some of Feattures from your EEG-Data during the Analysis of (STFT, Filtering and Windowing), so maybe you can try the HILBERT-HUANG-TRNASFORM and the EMPIRICAL Method to ectract the Frequancies-Bands. Look at this Link : http://www.mathworks.com/matlabcentral/fileexchange/19681-hilbert-huang-transform
Shitanshu Mishra
2012-3-20
Hello Kosai. Please help me in implementing your code. I am not able to implement your code. Can you please explain the code above. It would be great help. Please give atleast a little explaination about how to use your code on EEG RAW data. Its really very important for me! :/
Kosai
2012-3-20
hello Shitanshu Mishra, i will assume that you have the wavelet-toolbox in your Matlab version, because my code need it.
i will try to explain the code based on a (sin-signal);
1- we have to determine the Sampling Frequancy(Fs) in Hz:
FS = 1000;
2- Generate the Time Domain(t):
t = 1/Fs;
3- generate the SIN-Signal based on (t,Fs):
mySignal = sin(2*pi*t);
4- now we assume that mySignal relates to a EEG-RAW, and now we are going to do wavelet-analysis on this eeg-raw. as i mentioned before, the most important thing by wavelet-analysis is to determine the (Fs), in our case (Fs = 1000 hz). As you now the Bands Frequency for
Theta = 0-4 Hz,
Delta = 4-8 Hz,
Alpha = 8-16 Hz,
Beta = 16-32 Hz,
Gamma = 32-64 Hz,
and now to achive these Bands-Frequency, a wavelet-analysis is a very effective method for that, so wavelet analysis try to sink the frequency of our Signal(mySignal) in different Bands, so for every Bands we have to extract a Decompostion Level from the main Signal and do this continuously until we achive the Bands between (0-64 Hz), and now we have (Fs = 1000) then we need to decompose the Signal in (8) Levels to get the Bands between (0-64 Hz).
5- Determine the Wavelet Function to make the Analysis, there are to many functions-Family of wavelet, so i will advice you to use the (sym8 or the db8), in this example, i am going to use (db8):
waveletFunction = 'db8';
6- applicate the wavelet-analysis on mySignal :
[C,L] = wavedec(mySignal,8,waveletFunction);
7-Calculation The Coificients Vectors of every Band :
cD1 = detcoef(C,L,1); %NOISY
cD2 = detcoef(C,L,2); %NOISY
cD3 = detcoef(C,L,3); %NOISY
cD4 = detcoef(C,L,4); %NOISY
cD5 = detcoef(C,L,5); %GAMA
cD6 = detcoef(C,L,6); %BETA
cD7 = detcoef(C,L,7); %ALPHA
cD8 = detcoef(C,L,8); %THETA
cA8 = appcoef(C,L,waveletFunction,8); %DELTA
8- Calculation the Details Vectors of every Band :
D1 = wrcoef('d',C,L,waveletFunction,1); %NOISY
D2 = wrcoef('d',C,L,waveletFunction,2); %NOISY
D3 = wrcoef('d',C,L,waveletFunction,3); %NOISY
D4 = wrcoef('d',C,L,waveletFunction,4); %NOISY
D5 = wrcoef('d',C,L,waveletFunction,5); %GAMMA
D6 = wrcoef('d',C,L,waveletFunction,6); %BETA
D7 = wrcoef('d',C,L,waveletFunction,7); %ALPHA
D8 = wrcoef('d',C,L,waveletFunction,8); %THETA
A8 = wrcoef('a',C,L,waveletFunction,8); %DELTA
-------------
Shitanshu Mishra
2012-3-21
Thanks Kosai for such a favour.... My EEG Data had 256 Hz Sampling frequency, And I want to extract signals of all bands from given main signal. Please tell me that in how many levels should I decompose my Signal for Fs=256? Also tell me that how/where my final all 5 bands signals are returned, how can they be retrieved seperately? (I mean in what variable or function)?
Thanks for your help till now, Please do a little more favour. Regards. Shitanshu Mishra
Shitanshu Mishra
2012-3-21
Also:
Delta up to 4
Theta 4 – 8
Alpha 8 – 13
Beta >13 – 30
Gamma 30 – 10
Reference: http://en.wikipedia.org/wiki/Electroencephalography
Pleas help me how can I extract all 5 five bands Signals. I am new to Wavelet Decomposition.
dang huyen
2018-3-8
Hi Kosai, if I'm right, you are confusing here: detcoef syntax for detail coefficients, and wrcoef syntax for rescontruction. Should be reviewed again!
sunidhi gupta
2020-6-9
Hi kosai I tried ur code but it's showing error as " Expected so(data name) to be a vector" Can you please tell me how to rectify this error. Thank you in advance.
vinay kulkarni
2020-7-17
编辑：vinay kulkarni
2020-7-17
@sunidhi gupta I was facing same error. I added following line:
Data_Name=Data_Name(:);
And it works in my code.
It may help you.
Shitanshu Mishra
2012-3-21
0 个投票
Thanks Kosai for such a favour.... My EEG Data had 256 Hz Sampling frequency, And I want to extract signals of all bands from given main signal. Please tell me that in how many levels should I decompose my Signal for Fs=256? Also tell me that how/where my final all 5 bands signals are returned, how can they be retrieved seperately? (I mean in what variable or function)?
Thanks for your help till now, Please do a little more favour. Regards. Shitanshu Mishra
Shitanshu Mishra
2012-3-21
0 个投票
Also:
Delta up to 4
Theta 4 – 8
Alpha 8 – 13
Beta >13 – 30
Gamma 30 – 10
Pleas help me how can I extract all 5 five bands Signals. I am new to Wavelet Decomposition.
5 个评论
Kosai
2012-3-21
well, you have a EEG-Signal with Fs = 256 Hz , so in this case you need (5 Decomposition Levels) then :
function [Gamma,Beta,Alpha,Theta,Delta] = getBand(mySignal,waveletFunction)
[C,L] = wavedec(mySignal,5,waveletFunction);
%Calculation The Coificients Vectors of every Band
cD1 = detcoef(C,L,1); %NOISY
cD2 = detcoef(C,L,2); %Gamma
cD3 = detcoef(C,L,3); %Beta
cD4 = detcoef(C,L,4); %Alpha
cD5 = detcoef(C,L,5); %Theta
cA5 = appcoef(C,L,waveletFunction,5); %Delta
%Calculation the Details Vectors of every Band :
D1 = wrcoef('d',C,L,waveletFunction,1); %NOISY
D2 = wrcoef('d',C,L,waveletFunction,2); %Gamma
D3 = wrcoef('d',C,L,waveletFunction,3); %Beta
D4 = wrcoef('d',C,L,waveletFunction,4); %Alpha
D5 = wrcoef('d',C,L,waveletFunction,5); %Theta
A5 = wrcoef('a',C,L,waveletFunction,5); %Delta
Gamma = D2; figure, plot(1:1:length(Gamma),Gamma);
Beta = D3; figure, plot(1:1:length(Beta), Beta);
Alpha = D4; figure, plot(1:1:length(Alpha),Alpha);
Theta = D5; figure, plot(1:1:length(Theta),Theta);
Delta = A5; figure, plot(1:1:length(Delta),Delta);
% Call The Function :
% Fs = 256; t = 0:1/Fs:1; mySignal = sin(2*pi*t);
% waveletFunction = 'db8';
% [Gamma,Beta,Alpha,Theta,Delta]=getBand(mySignal,waveletFunction);
Kosai
2012-3-21
well, you have a EEG-Signal with Fs = 256 Hz , so in this case you need (5 Decomposition Levels) then :
function [Gamma,Beta,Alpha,Theta,Delta] = getBand(mySignal,waveletFunction)
[C,L] = wavedec(mySignal,5,waveletFunction);
%Calculation The Coificients Vectors of every Band
cD1 = detcoef(C,L,1); %NOISY
cD2 = detcoef(C,L,2); %Gamma
cD3 = detcoef(C,L,3); %Beta
cD4 = detcoef(C,L,4); %Alpha
cD5 = detcoef(C,L,5); %Theta
cA5 = appcoef(C,L,waveletFunction,5); %Delta
%Calculation the Details Vectors of every Band :
D1 = wrcoef('d',C,L,waveletFunction,1); %NOISY
D2 = wrcoef('d',C,L,waveletFunction,2); %Gamma
D3 = wrcoef('d',C,L,waveletFunction,3); %Beta
D4 = wrcoef('d',C,L,waveletFunction,4); %Alpha
D5 = wrcoef('d',C,L,waveletFunction,5); %Theta
A5 = wrcoef('a',C,L,waveletFunction,5); %Delta
Gamma = D2; figure, plot(1:1:length(Gamma),Gamma);
Beta = D3; figure, plot(1:1:length(Beta), Beta);
Alpha = D4; figure, plot(1:1:length(Alpha),Alpha);
Theta = D5; figure, plot(1:1:length(Theta),Theta);
Delta = A5; figure, plot(1:1:length(Delta),Delta);
% Call The Function :
% Fs = 256; t = 0:1/Fs:1; mySignal = sin(2*pi*t);
% waveletFunction = 'db8';
% [Gamma,Beta,Alpha,Theta,Delta]=getBand(mySignal,waveletFunction);
Shitanshu Mishra
2012-3-21
Hi Kosai!
You are so good. Thanks very much. You solved a big problem for me......
I am heartily thanking you. Thanks And warm wishes.
Shitanshu Mishra
Ari
2012-11-8
Hi Kosai, Thank you very much for your detailed post. I have a question regarding the frequency band when you use wavelet. The bands Shitanshu wanted are:
Delta up to 4 Theta 4 – 8 Alpha 8 – 13 Beta >13 – 30 Gamma 30 – 10
When we use wavedec, the wavelet of the signal is half at each level. So,we will have power computed for 8-16Hz, 16-32Hz, 32-64Hz instead of 8 – 13, 13 – 30Hz
Am I correct?
Also, for Fs=256, why do we have 5 decomposition levels. According to http://www.mathworks.com/help/wavelet/ref/wavedec.html
I think we should have 6 levels
Piyush kant
2017-10-20
Thanks for great explanation. I have a doubt, as you said Delta is 0-4Hz but in documentation it is given as 0.5-4Hz. the question is about that lower frequency, wouldn't that 0.0-0.5 Hz signal effect the remaining part of the same?
Salaheldin
2012-5-10
0 个投票
Hi Kosai,
Thank you very much for this very helpful code. I have two questions: (1) why do we need the cD=detcoef() and cA=appcoef() commands? We don't use cD or cA variables in any operations. (2) How does one calculate the power in the each of the decomposition levels, more specifically, bands of interest (gamma, beta, alpha, delta, and theta)?
one last point, if my sampling frequency is 2000, then I would need 9 decomposition levels to capture the bands of interest in the above expample, no?
thank you very much!!!
1 个评论
Kosai
2012-6-7
Hi Salaheldin,
actually we dont need the (cD,...,cA) but i have included them in code to have a better vision about Coificients and Details vectors in wavelet.
About the Power , simply calculate the Energy of every band (Energy = Power = abs(sum(myBand.^2)))
If your Fs is 2000 yes you need 9 Decomposition levels .
Hi KOSAI, Great explanation.I have a question. Suppose, my data is sampled with 500Hz. So, according to Nyqusit theorem there is 250 hz spectrum. I filtered the signal 0-64 Hz. I need 0-8Hz band. Now how much decomposition level is required 3 or 5?? Sampling freuency always 500 Hz.
Shishi
2014-3-27
0 个投票
Thank you Kosai, your piece of code and explanation was really helpful to me. As Riheen mentioned I'd like to add that the sampling frequency (Fs) is twice frequency spectrum based on Nyquist's theorem. Accordingly, if the sample frequency is 1000Hz then the frequency spectra is 500Hz and we need 7 levels of decomposition to get the delta band. You also can refer to the paper in below wherein the authors use the same approach. There is a figure of discrete wavelet packet transform (DWPT)which really shows the decomposition levels.
The sampling frequency of their EEG is 128 Hz.
VENKATA PHANIKRISHNA B
2017-5-15
编辑：VENKATA PHANIKRISHNA B
2017-5-15
decomposition level decision depends
on smaple frequency. it is ok. after decomposition
[C,L] = wavedec(mySignal,8,waveletFunction);
cD8 = detcoef(C,L,8); %THETA
cA8 = appcoef(C,L,waveletFunction,8); %DELTA
D8 = wrcoef('d',C,L,waveletFunction,8); %THETA
A8 = wrcoef('a',C,L,waveletFunction,8); %DELTA
among two which one is correct for getting Theta and Delta sub-bands of EEG signal. To calculate power and energy of theta and gamma sub-bands which signal is correct 'CD' or 'A'. Please provide the clarity.
Thank you
Jae Moon
2017-12-1
0 个投票
I've got an EEG data set that used a sampling freq of 512. I am trying to extract only the beta band (16-32Hz). What level of decomposition do I need?
1 个评论
Waqas Rasheed
2021-3-9
You may need 6, according to Kosai's logic.
Kosai! I tried to write code to find levels of decomposition following your logic. Please correct me if I'm wrong.
fs = 512; L = wmaxlev(fs,'db8'); LEVELS = L+round(L/4)
Gives 6 for 512, 8 for 1000, 9 for 2000, 5 for 256. Interestingly, it gives 4 for 128, and 3 for 64.
Now I have a question for Kosai, or anybody who can answer this. I have a signal at 5000 Fs, and from the above logic, it requires 10 levels of decomposition (if I am not wrong). Please let me know if teh following source is correct to extract delta, theta, alpha, beta, and gamma waves from any signal (say inputsignal).
waveletFunction = 'db8';
[C,L] = wavedec(inputsignal,10,waveletFunction);
GAMMA = wrcoef('d',C,L,waveletFunction,7);
BETA = wrcoef('d',C,L,waveletFunction,8);
ALPHA = wrcoef('d',C,L,waveletFunction,9);
THETA = wrcoef('d',C,L,waveletFunction,10);
DELTA = wrcoef('a',C,L,waveletFunction,10);
... and 1 - 6 are NOISE, right?
类别
在 帮助中心 和 File Exchange 中查找有关 Hilbert and Walsh-Hadamard Transforms 的更多信息
2021-3-9
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
