[Eeglablist] MATLAB fft vs spectopo output

[Makoto Miyakoshi]

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>