[Eeglablist] MATLAB fft vs spectopo output
Samuel Tomlinson
stomlin1 at swarthmore.edu
Wed Jun 25 10:22:58 PDT 2014
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>
More information about the eeglablist
mailing list