[Eeglablist] MATLAB fft vs spectopo output
Makoto Miyakoshi
mmiyakoshi at ucsd.edu
Fri Jun 27 21:36:02 PDT 2014
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20140627/be343100/attachment.html>
More information about the eeglablist
mailing list