[Eeglablist] exact same values of spectopo.m using pwelch()

Makoto Miyakoshi mmiyakoshi at ucsd.edu
Thu Dec 25 19:01:57 PST 2014


Dear James,

Without trying your code, let me tell you that Andreas Widmann once
clarified this issue on this list. Let me guide your attention to it--see
below. You should also be able to find it in the eeglablist archive.

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

On Thu, Dec 18, 2014 at 4:32 PM, James Bonello <jamesbonello9 at gmail.com>
wrote:

> Hi EEGLablist,
>
> The spectopo function calculates power spectral density values for a given
> signal. However I would like to try get the same exact values of the
> spectopo function but using MATLAB's pwelch() function instead since this
> function should give the power spectral density as well. However I am
> getting completely different values and I am not sure if the way I am
> calculating it is correct:
>
> X  = EEG.data(channel_no,:,epoch_no);
>          length_X  = EEG.pnts;                    % signal length
>          nfft = 2^nextpow2(length_X);             % Next power of 2 from
> length of x
>
>          no_of_overlapped_samples = 0;
>          window = 512;
>
>          % get psd using pwelch()
>          [spectra,freqs] = pwelch(X, window, no_of_overlapped_samples,
> nfft, EEG.srate);
>
>          spectra = 10*log(spectra);
>
> Above I am calculating it using the pwelch() function, with spectopo it is
> as below:
>
> [spectra,freqs] = spectopo(EEG.data(channel_no,:,epoch_no), 0, EEG.srate);
>
>
> Is there something I am missing or am doing wrong?
>
> Thank you!
>
> --
> *James Bonello*
>
> _______________________________________________
> 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/20141225/137ca943/attachment.html>


More information about the eeglablist mailing list