[Eeglablist] MATLAB fft vs spectopo output
Andreas Widmann
widmann at uni-leipzig.de
Sat Jun 28 07:18:52 PDT 2014
Hi Makoto and Sam,
on the risk also getting lost in the energy/power/spectrum/spectral density-jungle:
* Units were mixed: spectopo/pwelch compute PSD, that is µV^2/Hz. So you have to normalize by sampling rate not data/fft length.
* Data were windowed but not normalized to window power.
* The whole vector was divided by 2 two times instead of only the first and last point.
Find some sample code how two estimate PSD from FFT here:
http://www.mathworks.com/matlabcentral/answers/33653-psd-estimation-fft-vs-welch
If you use the code below you should get perfectly identical results (the modified lines are marked by capital comments).
Hope this helps! Best,
Andreas
% do it in individual steps for debugging
w = hamming(N);
singleTrial = singleTrial.*w';
pow = fft(singleTrial); % apply fft
pow = pow(1:floor(N/2+1)); % toss redundant half of dataset
pow = abs(pow); % absolute value for magnitude
pow = pow.^2; % square to get power
pow = pow/(w'*w); % NORMALIZE BY WINDOW POWER
pow = pow./srate; % NORMALIZE BY SAMPLING RATE
% pow = pow./N; % normalize by length of data
pow = pow.*2; % mult power by 2 b/c removed half of dataset
% DC should be unique-- undo mult. by 2
% pow(1,:) = pow(1,:)./2;
size(pow) % !!!
pow(1) = pow(1)./2; % ONLY FIRST AND LAST POINT (NOT THE VECTOR)
% so should nyquist-- undo mult. by 2 if nyquist is included
if ~rem(N,2)
% here nyquist is included- should be unique as well
% pow(end,:) = pow(end,:)./2;
pow(end) = pow(end)./2; % ONLY FIRST AND LAST POINT (NOT THE VECTOR)
end
Am 28.06.2014 um 06:36 schrieb Makoto Miyakoshi <mmiyakoshi at ucsd.edu>:
> Dear Samuel,
>
> Hmm interesting. Probably it's an issue of window function difference between Hamming and the one used in welch which I could not confirm... Andreas, why don't you try Samuel's code, it works by copy and paste on Matlab. I would appreciate your comment.
>
> -1 to 1 uV peak-to-peak sine wave measures 0dB at EEGLAB spectopo(), which makes sense to me.
>
> Maktoo
>
>
>
>
> On Wed, Jun 25, 2014 at 10:22 AM, Samuel Tomlinson <stomlin1 at swarthmore.edu> wrote:
> Hello eeglablist,
>
> I am trying to reconcile the results from calculating an FFT from scratch with the results I get using EEGLAB's spectopo function. I'm at the point where I'm estimating absolute power from an epoched dataset with MATLAB's fft and applying a hamming window to match the pwelch function (which is called by spectopo).
>
> I've attached some code below that compares the results of this procedure with the output from spectopo. The code generates a plot of the power spectrum in dB. As you'll see, the rank correlation between the fft method and the spectopo method is very high (~0.99). But for some reason, the fft results are shifted down (linearly, it seems) by 15 units compared to the spectopo results. I've worked through the code for spectopo and pwelch and can not account for this apparent shift.
>
> I know this is a vague question, but has anyone else experienced scaling issues between spectopo's output and the results from built-in MATLAB functions? Any guidance would be much appreciated.
>
> Thank you,
> Sam Tomlinson
>
>
> Here's the code:
>
> % Read in data- nchans x npoints x nepochs for one experimental condition
> data = rand(128,750,20); % run on random data for demonstration
> srate = 250; % 250 Hz, 3 second epochs, 20 epochs
>
> % get some params from EEG dataset
> nchans = size(data,1);
> N = size(data,2);
> nepochs = size(data,3);
>
> % set some standard values for fft
> nyquist = srate/2;
> freqs = linspace(0,nyquist,floor(N/2+1));
>
> % initialize output matrix
> eegspec = zeros(nchans,floor(N/2+1));
>
> % loop thorugh channels and compute spectral power- update eegspec
> for chan = 1:nchans
>
> % matrix to store fft output for each epoch
> epochData = zeros(floor(N/2+1),nepochs);
>
> % loop epochs for each channel
> for epoch = 1:nepochs
>
> singleTrial = data(chan,:,epoch);
> % DC offset
> % singleTrial = singleTrial - mean(singleTrial);
>
> % apply window, perform fft
> % pow = abs(fft(x.*hann(N)')/N).^2;
>
> % do it in individual steps for debugging
> singleTrial = singleTrial.*hamming(N)';
> pow = fft(singleTrial); % apply fft
> pow = pow(1:floor(N/2+1)); % toss redundant half of dataset
> pow = abs(pow); % absolute value for magnitude
> pow = pow./N; % normalize by length of data
> pow = pow.^2; % square to get power
> pow = pow.*2; % mult power by 2 b/c removed half of dataset
>
> % DC should be unique-- undo mult. by 2
> pow(1,:) = pow(1,:)./2;
> % so should nyquist-- undo mult. by 2 if nyquist is included
> if ~rem(N,2)
> % here nyquist is included- should be unique as well
> pow(end,:) = pow(end,:)./2;
> end
>
> epochData(:,epoch) = pow;
>
> end
>
> epochData = mean(epochData,2); % mean of rows-- avg across epochs
> eegspec(chan,:) = epochData';
>
> end
>
> % % way to condense lines 19-55
> % 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
>
> freqs = freqs(1:size(eegspec,2)); % make sure dimensions match
>
>
> % ######## run comparison with EEGLAB spectopo ######
>
> chan = randi([1,128]); % get a random channel
>
> % force spectopo into same parameters that we use for fft; specify zero padding
> [a,b] = spectopo(data,750,250,'freqrange',[1 50],'nfft',750,'winsize',750,'freqfac',1);
> close(gcf)
>
> % first, show that frequencies are same
> disp('frequency: ')
> corrcoef(freqs,b)
>
> % now, compare raw (unlogged) power (microV^2)
> a_unlog = 10.^(a/10);
> disp('spectrum unlogged: ')
> corrcoef(eegspec(chan,:),a_unlog(chan,:))
>
> %eegspec = eegspec.*sqrt(N);
>
> % now, compare logged (dB) power spectra
> eegspec_log = 10.*log10(eegspec);
> disp('spectrum dB (10*log10): ')
> corrcoef(eegspec_log(chan,:),a(chan,:))
>
> % now, find 50 Hz index for plotting, plot our fft results with spectopo
> % for the random channel
> rowsInx = find(freqs<=50);
> upperLim = max(rowsInx);
>
> plot(freqs(1:upperLim),eegspec_log(chan,1:upperLim),'LineWidth',2,'Color','r'); % ours
> hold on
> plot(freqs(1:upperLim),a(chan,1:upperLim),'LineWidth',2,'Color','b') % spectopo's
> ylabel(gca,'Power dB == 10*log_{10}(\muV^{2})');
> xlabel(gca,'Frequency (Hz)');
>
>
>
> _______________________________________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.ucsd.edu
> For digest mode, send an email with the subject "set digest mime" to eeglablist-request at sccn.ucsd.edu
>
>
>
> --
> Makoto Miyakoshi
> Swartz Center for Computational Neuroscience
> Institute for Neural Computation, University of California San Diego
> _______________________________________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.ucsd.edu
> For digest mode, send an email with the subject "set digest mime" to eeglablist-request at sccn.ucsd.edu
More information about the eeglablist
mailing list