[Eeglablist] MATLAB fft vs spectopo output
Makoto Miyakoshi
mmiyakoshi at ucsd.edu
Mon Jun 30 09:54:49 PDT 2014
Quick and accurate, AWESOME!
Thank you Andreas.
Makoto
On Sat, Jun 28, 2014 at 7:18 AM, Andreas Widmann <widmann at uni-leipzig.de>
wrote:
> 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
>
>
--
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/20140630/61bb3123/attachment.html>
More information about the eeglablist
mailing list