FFT vs eeglab spectopo

[adamn]

Hi everyone, 

I'm trying to reconcile the results from calculating an FFT as shown in the book, and the outputs I get from EEGLAB using the spectopo function. 

I'm at the point where I'm using a single channel in a single trial, and I'm having spectopo return the same number of frequency bins (no padding). The plots of the resulting spectra are similar, but the data aren't the same (rank correlation of ~ .7).

I'm estimating power using the looping with sine waves and the FFT command and getting the same results. I'm also using this approach for scaling 

http://sccn.ucsd.edu/pipermail/eeglablist/2010/003550.html

which yields different absolute numbers, but they correlate perfectly with the FFT results.

But I still don't get the same results. 

As far as I can tell, there are 2 things happening here that I don't yet fully grasp. 

(1) EEGLAB is doing some transformation of the values

(2) EEGLAB is using pwelch, not FFT, which seems to take a different approach for calculating power. 

I realize that this is a vague question, but I'm hoping that someone can shed some light on this inconsistency in different methods of estimating power. 

Also, the book is really really fantastic.

Thanks!

-Adam

[Mike X Cohen]

Hi Adam. 

I'm not overly familiar with eeglab's spectopo function, so off the top of my head, I don't know what would be producing the difference. However, a correlation of 0.7 sounds like nothing to worry about to me. There are several parameters to a short-time FFT approach, and any minor deviations will produce different results. Even something minor like changing the percent overlap of the time windows will already change the results slightly. 

In general, I would say -- purely based on your correlation of 0.7 -- that this is nothing to worry about. If you would like to understand what exactly is causing this difference, from a didactic perspective, then I'll need some more information. Perhaps if you send me your script (you could email it to me personally rather than posting it on the list, if you prefer) I can dig into it a bit. 

Mike

--
You received this message because you are subscribed to the Google Groups "AnalyzingNeuralTimeSeriesData" group.
To unsubscribe from this group and stop receiving emails from it, send an email to analyzingneuraltimes...@googlegroups.com.
Visit this group at http://groups.google.com/group/analyzingneuraltimeseriesdata.
For more options, visit https://groups.google.com/d/optout.

--

Mike X Cohen, PhD
mikexcohen.com

[Mike X Cohen]

Hi Adam,

(To update the list: Adam sent me some data and a script off-list.) There are four differences between what you did [2*abs(fft(data)/n)] and what the eeglab function spectopo does. 

1) spectopo, or more specifically, computeperiodogram.m, which is called by welch.m, which is called by pwelch.m, which is called by spectopo, is computing power (abs(fft(data)/n).^2, which is actually computed as the equivalent-but-slightly-faster complex conjugate (line 73). See page 345 in the book. You computed amplitude, which is the square root of power. 

2) spectopo.m computes 'raw' power but returns the dB power (see line 343: eegspecdB = 10*log10(eegspecdB); )

3) This one is weird: spectopo.m is actually only using the first 256 data points and ignoring the rest (it's still computing a 500-point FFT, so there is zero-padding being done). See lines 861-864: winlength ends up being 256. I'm not sure why eeglab has it like this. 

4) spectopo.m is tapering the data using a Hamming window (actually performed in welchparse.m on line 160). Your script doesn't apply a taper. 

So there you have it. Basically, you are computing amplitude and eeglab is returning dB power from the tapered first half of your data. 

Personally, I would recommend using your code instead of eeglab's. I would also recommend computing power and tapering the data. You can use the following code: 

pow = abs( fft(data.*hann(length(data))')/N ).^2;

I hope that's clear. 

Mike

On Wed, Jun 18, 2014 at 7:29 PM, Naples, Adam <adam....@yale.edu> wrote:

Thank you so much for helping us with this. 

I'm also copying a student in the lab on this email who is working with me on this.

What concerned me was most, was that this was not a short time FFT, but just an FFT (we're trying to look at resting alpha asymmetry). 

I figured that we might see some differences in rounding, but that overall, the results would correlate more than .7. 

I'm attaching my sample data and my code (The data is a 128x500*18 matrix called EEG). 

My hope is that what is happening is that I'm scaling the results of the FFT in some nonlinear and inappropriate way, but I'm not sure if my way or the EEGLAB way is appropriate. I haven't been able to get a good grasp of the different ways to transform FFT results to estimate power, and I haven't found a good resource for understanding this step. 

I also wonder if eeglab is doing some smoothing to the data? aside the the differences in scale, the plot of the spectral profile from eeglab is much smoother.

I'm attaching my code and the MAT file. 

Thanks again very much. We've already bought two copies of your book in the lab. It is tremendously clear and helpful. 

Best, 

Adam

[adamn]

Hi again, we're still trying to sort out our spectral power approach using EEGlab's spectopo as a point of reference. I think that our current concern is that we're scaling the results of our FFT incorrectly, some of our values seem very low compared to what we've seen in the literature and the results of EEGLAB, however the results correlate .99 with EEGLAB. 

This leads me to believe that we're missing something when scaling the data out of our FFT. 

Once again, I'm attaching our code. Most of the code is looping through channels and epochs and comparing it against the EEGLAB results. The meat of the code is in lines 32-38. 

Any guidance would be great. we're trying to calculate ratios between different power bands, but when one of the bands has power almost at 0 our results fluctuate quite a bit. 

Thanks

-a

[Mike X Cohen]

Hi Adam. The code looks fine to me. It crashed with the correlations because the frequency vectors were not the same. The values also seem reasonable to me -- alpha power of around 2-5 units, which I assume are microvolts. The surrounding frequencies do indeed have low power, but is perhaps due to having a small amount of clean data with a strong alpha rhythmicity. 

btw, your lines 19-55 could be condensed to the following: 

% loop thorugh channels and compute spectral power- update eegspec

for chan = 1:nchans

    pow = mean(abs(fft(bsxfun(@times,squeeze(data(chan,:,:)),hann(N)))/N).^2,2);

    eegspec(chan,:) = pow(1:length(freqs));

end

In short, I think your code is OK. Perhaps try running it on more data (I assume you have more than 500 time points). 

Mike

--
You received this message because you are subscribed to the Google Groups "AnalyzingNeuralTimeSeriesData" group.
To unsubscribe from this group and stop receiving emails from it, send an email to analyzingneuraltimes...@googlegroups.com.
Visit this group at http://groups.google.com/group/analyzingneuraltimeseriesdata.
For more options, visit https://groups.google.com/d/optout.

[Matt Euler]

Dear Dr. Cohen,

I have found this thread very helpful for trying to troubleshoot essentially the same question that Adam had. I have just a quick question about the formula for obtaining power that you provided:

pow = abs( fft(data.*hann(length(data))')/N ).^2;

Since we don't need to run the inverse fft at a later point in this application, isn't it also necessary to multiply by 2 at some point (excluding the DC) to preserve the signal power that would otherwise be contained in the negative frequencies? (e.g., per your FunTF book, section 4.3). I realize that a linear operation like that won't affect patterns of results, just their scale, but I'm just wanting to be sure I understand the method as fully as possible. 

Thanks,

Matt Euler

To unsubscribe from this group and stop receiving emails from it, send an email to analyzingneuraltimeseriesdata+unsub...@googlegroups.com.

Visit this group at http://groups.google.com/group/analyzingneuraltimeseriesdata.
For more options, visit https://groups.google.com/d/optout.

[Mike X Cohen]

Hi Matt. Indeed, the neuroscience community tends to be a bit fast-and-loose with units. In part I think that's OK because absolute units are difficult to interpret and even more difficult to compare across studies and methods. You are correct that there should be a *2 in that formula for an accurate reconstruction. And also as you mention, it doesn't affect the overall pattern (although squaring vs. not squaring does affect the shape of the power spectrum).

I discuss this issue in more detail in the lecture called "Accurately recovering data units in FFT and convolution" on my website.

Mike

Visit this group at https://groups.google.com/group/analyzingneuraltimeseriesdata.

For more options, visit https://groups.google.com/d/optout.

[Courtney J]

Hi all, 

 

We have a similar problem to those above in regard to reconciling results from pwelch and spectopo. We have tried ensuring our parameters are the same, and we have converted spectopo into mv. 

 

Our spectopo code is: 

 

[spectra,freqs] = spectopo(EEG.data,length(EEG.data),EEG.srate,'winsize',2500,'nfft',500,'overlap',1250); % Match Spectopo Parameters to PWELCH

spectra_mv = 10.^(spectra/10); % Convert Back to mV from dB

spectra_dmv = log10(1+spectra_mv); % Log Transform the Resulting PSD Matrix to Obtain Better Distribution

spectra_tdmv = spectra_dmv'; % Transpose Data

spectra_ftdmv = spectra_tdmv(1:40,:); % Extract FREQRANGE 1:40

clear spectra freqs spectra_mv spectra_dmv spectra_tdmv

 

Our personal pwelch script uses:

EEG_test1093 = EEG.data'; % transpose 

restEEG = EEG_test1093; % identify the experimental data you are working with

fs=500; %sampling rate

window=2500; %at a sampling rate of 500hz, 2500samples is 5 sec

freqrange=1:40; %identify the frequency range you want to extract PSD in

 

computing the ERD across 1:fqs

% calculate PSDs useing pwelch.m

[restPSD,freqs]=pwelch(restEEG,window,[],freqrange,fs); % the [] is to signal the default of 50% overlapping window

 

% log transform the resulting PSD matrix to obtain better distribution

restPSDlog=log10(1+restPSD); % takes the natural log

PSDlog_test1093 = restPSDlog;

clear fs EEG_1093RS restEEG freqrange restPSD restPSDlog window

 

When searching in spectopo, we noticed that pwelch is used in a for-loop to perform multiple pwelch calculations (spectopo lines 912-974)  - however our script performs only one calculation. Additionally, it seems as though spectopo has deconstructed the EEG.data structure, and reconstructs it after performing the pwelch calculations (using the matsel function in the pwelch line).  Lastly, we have noticed that transposing data occurs at different times in the calculations. 

 

Most importantly, is our pwelch calculation correct? And if possible, can you help us understand why the results are different? 

 

We have attached the code.

 

Much thanks in advance. 

 

Courtney 

[Mike X Cohen]

Hi Courtney. I don't immediately see anything wrong with your code, but sometimes it's difficult to tell without data. Perhaps you can share some screenshots of what the plots look like from spectopo vs. pwelch? Looking back through this email thread shows a message from June 2014 in which I identified several differences between spectopo and fft. I'm not sure if spectopo has been adjusted since then.

Mike

Mike X Cohen, PhD
New online courses: mikexcohen.com

[Courtney J]

Hi Mike, 

Thank you for your prompt response!

I have attached plots for both output variables. 

We have attempted to account for the differences explained above by: 

1) using pwelch instead of fft in the pwelch code

2) converting the spectopo output back to mV

3) specifying winlength (winsize) in the spectopo code 

4) and using the default Hamming window in the pwelch code

The only other difference we found was that spectopo uses a for-loop (for nchans) to perform multiple pwelch calculations on a data structure with only one channel  - however our script performs only one calculation on the entire EEG.data structure (nchans, frames, epochs).

Any additional help would be much appreciated in confirming our power calculation results or reconciling the spectopo difference. Thanks again very much for your time!

Courtney

[Mike X Cohen]

Hi Courtney. If you set those two figures to the same y-axis limits, you'll see that they are extremely really super-duper a lot similar to each other. It seems like the main difference is a bit of 1/f scaling. If you are interested in figuring out the exact cause of the difference as an exercise in MATLAB programming, then go for it. Otherwise, I would consider the two results to be the same. 

Mike

[Steve Miller]

PWELCH looks like magnitude while SPECTOPO looks like power, dunno if that's what you meant.