[Eeglablist] MATLAB fft vs spectopo output

[Samuel Tomlinson]

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)'); 


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20140625/24582d83/attachment.html>