Hi Diogo,

are you using "Channel spectra and maps" to get your plots? If so, you will want to use the underlying spectopo-function to get your mean values. In your case, the commands may look something like this (assuming theta = 4-8 Hz):

sampRate = xxx; %sampling rate of your data

lowerFreq = 4; %lower bound of your frequency band of interest

higherFreq = 8; %upper bound of your frequency band

chanNr = xx; %channel number of your specific lead

%computing log spectrum for different frequencies

[power, freq] = spectopo(EEG.data(chanNr, :), 0, sampRate);

%average power within the predefined frequency range

meanTheta = mean(power(find(freq == lowerFreq) : find(freq == higherFreq)));

I hope, this will fit your needs. Otherwise, please let me know.

Best

Chris

Hi Cris,

Thank you very much for helping me.

I used this code but the variable "meantheta" in workspace resulted in NaN.

There is another variable,"power" 1x513. May I use these values to cumpute the mean?

Best,

Hi Diogo,

turned out that my script was not flexible enough to deal with certain sampling rates. I assume that was the problem. I send you the new version in which actually only the last line has changed. By the way: you are right about the power vector. However, the power vector contains power values of all computable frequencies. So you still need to extract a subset for your frequency band of interest. Basically, the freq-vector contains the labels for the x-axis, the power-vector contains the y-values. However, everything should go automatically now.

sampRate = xxx; %sampling rate of your data

lowerFreq = 4; %lower bound of your frequency band of interest

higherFreq = 8; %upper bound of your frequency band

chanNr = xx; %channel number of your lead

%computing log spectrum for different frequencies

[power, freq] = spectopo(EEG.data(chanNr, :), 0, sampRate);

%average power within the predefined frequency range

meanPower = mean(power(freq >= lowerFreq & freq

Hi Cris,

It´s working fine now!

Thank a lot once again!

Best

Diogo

Hi Federico,

thank you for sending me some sample data. Unfortunately, I did not run (automatically) into the same problem. The frequency resolution depends on the size of your FFT window. Normally, the spectopo-function automatically chooses adequate parameters for the window size based on your data. You could try to set the number of data sample points per window manually and look if this was the problem. A suggestion (from 'my' eeglab for your data) would be 512 samples, you just need to add the winsize argument in the second last line:

[power, freq] = spectopo(EEG.data(chanNr, :), 0, sampRate, 'winsize', 512);

As the spectopo-function oversamples (x2) the data prior to the FFT by default, you should get an actual FFT window size of 1024 samples and a frequency resolution of .4883 Hz / bin (500/1024). To be safe you can set the oversampling factor manually, too.

[power, freq] = spectopo(EEG.data(chanNr, :), 0, sampRate, 'winsize', 512, 'freqfac', 2);

Please let me know if this works out for you.

Best,

Chris

Hi Chris,

I had the very same doubt Diogo did and your script was a good solution for me. Thanks a lot for sharing! =)

However, I have a question: when you run this script does it take into account the whole epoch, right? I am asking this because I have a 500ms baseline period before the interest trigger that I would like to exclude from PSD calculation... Do I have to exclude mannually this baseline period in order to extract values I am interested?

Best,

Nastassja

P.s.: do you know how to extract ERSP values as well?

Hi Chris,

It is really great feedback and explanation from you.

I also have other question.

I have dataset with 14 channels of raw eeg data (like 14x1190).

I want to filter this one with low and high bandwith (1~43, delta~gamma) and get PSD values and plot that data on matlab. Could you give me some advice or example that I can visualize the data like scatter format? I am pretty new to this subject, so it is so difficult for me. I know how to plot them through eeglab, but not getting actual PSD data after filter out.

Thank you

Hey Nastassja,

sorry that I did not reply to you earlier. I somehow must have missed the notification of your post. I am not sure if it is still of use for you but here would be my answer:

You're right, the script I posted uses the whole available time window. One possibility would be if you cut off everything around. However, I think that is more convenient to just indicate the sample points you want to analyze. Below you'll find the previous script with an extra line where you can enter the sample points that form your time window of interest using the sampWin variable (the values below are just some arbitrary numbers).

sampRate = xxx; %sampling rate of your data

lowerFreq = 4; %lower bound of your frequency band of interest

higherFreq = 8; %upper bound of your frequency band

chanNr = xx; %channel number of your lead

sampWin = 101:500; %indices of sample points you want to go into analysis. If empty: whole segment is used.

%computing log spectrum for different frequencies

if isempty(sampWin)

[power, freq] = spectopo(EEG.data(chanNr, :), 0, sampRate);

else

[power, freq] = spectopo(EEG.data(chanNr, sampWin), 0, sampRate);

end

%average power within the predefined frequency range

meanPower = mean(power(freq >= lowerFreq & freq

Hi Yonggun,

I am not 100% sure if I fully understood your question so please tell me if my reply doesn't answer it. The getting-PSD-values part is basically what Diogo asked, so it should be possible to use the script posted above. About the scatter plot. part: I assume that you want to do a scatter plot with frequency on the x axis and power on the y axis? Then you just could use the power and freq vectors that the script returns. freq on the x axis, and power on the y axis:

figure, scatter(freq, power);

Hope that helps you. Good luck,

Chris

Do you can explain how to use it powerband, They are the features for the classification ? I'm working with ERS/ERD, Following Bernhard Graimann and Gert Pfurtscheller (attach file)  but I don't understand the feature extraction yet... Thanks. 

 Hi Chris,

I have one more question related to the discussion above. 

I'm trying to run the code you mention to calculate power values, but I'm getting the following error: Subscript indices must either be real positive integers or logicals.

Does that mean that I'm getting negative values for the power? 

One more question: Are the power values that EEGlab produce get corrected to baseline? and how do I tell EEGlab which is the baseline period I'm interested in?

Thank you 

Hi Chris,

I’m also new to PSD and was wondering how you could use this script on a 20 minute single channel EEG recording to analyze the change in alpha level activity over time, instead of just the mean?

Thanks,

Thomas

I could also use a little help. I'm on a similar vein, and have been able to use some lines of the code but not all. 

Ideally, I'd like to have average power head maps for all bands across my 19point data. I'm more familiar with the eeglab gui than matlab commands but can probably make sense of either.

Indeed this is a fantastic code! The output seems to be in dB. What would be the function to get an output in µV²?

Hi can you use Darbelai It is an extension in for EEGlab. You will find it under manage extensions. I have used this to get power for theta, alpha, etc. I attached some screen shots to get you to the extension. Maybe that will help you!

Hi, this is great code! I’m getting an output that is almost what I needed, but overlap is 0% and I need 10%. I’m using the code [power,freq] = spectopo(EEG.data, 0, 1000);. Where in that code would I add in the 10% overlap?

thanks!