<div dir="ltr">Dear James,<div><br></div><div>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.</div><div><br></div><div>Makoto</div><div><br></div><div><br class=""><br><div class="gmail_quote">On Sat, Jun 28, 2014 at 7:18 AM, Andreas Widmann <span dir="ltr"><<a href="mailto:widmann@uni-leipzig.de" target="_blank">widmann@uni-leipzig.de</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Hi Makoto and Sam,<br><br>on the risk also getting lost in the energy/power/spectrum/spectral density-jungle:<br>* Units were mixed: spectopo/pwelch compute PSD, that is µV^2/Hz. So you have to normalize by sampling rate not data/fft length.<br>* Data were windowed but not normalized to window power.<br>* The whole vector was divided by 2 two times instead of only the first and last point.<br><br>Find some sample code how two estimate PSD from FFT here:<br><a href="http://www.mathworks.com/matlabcentral/answers/33653-psd-estimation-fft-vs-welch" target="_blank">http://www.mathworks.com/matlabcentral/answers/33653-psd-estimation-fft-vs-welch</a><br><br>If you use the code below you should get perfectly identical results (the modified lines are marked by capital comments).<br><br>Hope this helps! Best,<br>Andreas<br><span class=""><br> % do it in individual steps for debugging<br></span> w = hamming(N);<br> singleTrial = singleTrial.*w';<br><span class=""> pow = fft(singleTrial); % apply fft<br> pow = pow(1:floor(N/2+1)); % toss redundant half of dataset<br> pow = abs(pow); % absolute value for magnitude<br></span><span class=""> pow = pow.^2; % square to get power<br></span> pow = pow/(w'*w); % NORMALIZE BY WINDOW POWER<br> pow = pow./srate; % NORMALIZE BY SAMPLING RATE<br>% pow = pow./N; % normalize by length of data<br><span class=""> pow = pow.*2; % mult power by 2 b/c removed half of dataset<br><br> % DC should be unique-- undo mult. by 2<br></span>% pow(1,:) = pow(1,:)./2;<br> size(pow) % !!!<br> pow(1) = pow(1)./2; % ONLY FIRST AND LAST POINT (NOT THE VECTOR)<br><span class=""> % so should nyquist-- undo mult. by 2 if nyquist is included<br> if ~rem(N,2)<br> % here nyquist is included- should be unique as well<br></span>% pow(end,:) = pow(end,:)./2;<br> pow(end) = pow(end)./2; % ONLY FIRST AND LAST POINT (NOT THE VECTOR)<br> end<br><br><br>Am 28.06.2014 um 06:36 schrieb Makoto Miyakoshi <<a href="mailto:mmiyakoshi@ucsd.edu">mmiyakoshi@ucsd.edu</a>>:<br><div class=""><div class="h5"><br>> Dear Samuel,<br>><br>> 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.<br>><br>> -1 to 1 uV peak-to-peak sine wave measures 0dB at EEGLAB spectopo(), which makes sense to me.<br>><br>> Maktoo<br>><br>><br>><br>><br>> On Wed, Jun 25, 2014 at 10:22 AM, Samuel Tomlinson <<a href="mailto:stomlin1@swarthmore.edu">stomlin1@swarthmore.edu</a>> wrote:<br>> Hello eeglablist,<br>><br>> 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).<br>><br>> 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.<br>><br>> 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.<br>><br>> Thank you,<br>> Sam Tomlinson<br>><br>><br>> Here's the code:<br>><br>> % Read in data- nchans x npoints x nepochs for one experimental condition<br>> data = rand(128,750,20); % run on random data for demonstration<br>> srate = 250; % 250 Hz, 3 second epochs, 20 epochs<br>><br>> % get some params from EEG dataset<br>> nchans = size(data,1);<br>> N = size(data,2);<br>> nepochs = size(data,3);<br>><br>> % set some standard values for fft<br>> nyquist = srate/2;<br>> freqs = linspace(0,nyquist,floor(N/2+1));<br>><br>> % initialize output matrix<br>> eegspec = zeros(nchans,floor(N/2+1));<br>><br>> % loop thorugh channels and compute spectral power- update eegspec<br>> for chan = 1:nchans<br>><br>> % matrix to store fft output for each epoch<br>> epochData = zeros(floor(N/2+1),nepochs);<br>><br>> % loop epochs for each channel<br>> for epoch = 1:nepochs<br>><br>> singleTrial = data(chan,:,epoch);<br>> % DC offset<br>> % singleTrial = singleTrial - mean(singleTrial);<br>><br>> % apply window, perform fft<br>> % pow = abs(fft(x.*hann(N)')/N).^2;<br>><br>> % do it in individual steps for debugging<br>> singleTrial = singleTrial.*hamming(N)';<br>> pow = fft(singleTrial); % apply fft<br>> pow = pow(1:floor(N/2+1)); % toss redundant half of dataset<br>> pow = abs(pow); % absolute value for magnitude<br>> pow = pow./N; % normalize by length of data<br>> pow = pow.^2; % square to get power<br>> pow = pow.*2; % mult power by 2 b/c removed half of dataset<br>><br>> % DC should be unique-- undo mult. by 2<br>> pow(1,:) = pow(1,:)./2;<br>> % so should nyquist-- undo mult. by 2 if nyquist is included<br>> if ~rem(N,2)<br>> % here nyquist is included- should be unique as well<br>> pow(end,:) = pow(end,:)./2;<br>> end<br>><br>> epochData(:,epoch) = pow;<br>><br>> end<br>><br>> epochData = mean(epochData,2); % mean of rows-- avg across epochs<br>> eegspec(chan,:) = epochData';<br>><br>> end<br>><br>> % % way to condense lines 19-55<br>> % for chan = 1:nchans<br>> % pow = mean(abs(fft(bsxfun(@times,squeeze(data(chan,:,:)),hann(N)))/N).^2,2);<br>> % eegspec(chan,:) = pow(1:length(freqs));<br>> % end<br>><br>> freqs = freqs(1:size(eegspec,2)); % make sure dimensions match<br>><br>><br>> % ######## run comparison with EEGLAB spectopo ######<br>><br>> chan = randi([1,128]); % get a random channel<br>><br>> % force spectopo into same parameters that we use for fft; specify zero padding<br>> [a,b] = spectopo(data,750,250,'freqrange',[1 50],'nfft',750,'winsize',750,'freqfac',1);<br>> close(gcf)<br>><br>> % first, show that frequencies are same<br>> disp('frequency: ')<br>> corrcoef(freqs,b)<br>><br>> % now, compare raw (unlogged) power (microV^2)<br>> a_unlog = 10.^(a/10);<br>> disp('spectrum unlogged: ')<br>> corrcoef(eegspec(chan,:),a_unlog(chan,:))<br>><br>> %eegspec = eegspec.*sqrt(N);<br>><br>> % now, compare logged (dB) power spectra<br>> eegspec_log = 10.*log10(eegspec);<br>> disp('spectrum dB (10*log10): ')<br>> corrcoef(eegspec_log(chan,:),a(chan,:))<br>><br>> % now, find 50 Hz index for plotting, plot our fft results with spectopo<br>> % for the random channel<br>> rowsInx = find(freqs<=50);<br>> upperLim = max(rowsInx);<br>><br>> plot(freqs(1:upperLim),eegspec_log(chan,1:upperLim),'LineWidth',2,'Color','r'); % ours<br>> hold on<br>> plot(freqs(1:upperLim),a(chan,1:upperLim),'LineWidth',2,'Color','b') % spectopo's<br>> ylabel(gca,'Power dB == 10*log_{10}(\muV^{2})');<br>> xlabel(gca,'Frequency (Hz)');<br>><br>><br>><br>> _______________________________________________<br>> Eeglablist page: <a href="http://sccn.ucsd.edu/eeglab/eeglabmail.html" target="_blank">http://sccn.ucsd.edu/eeglab/eeglabmail.html</a><br>> To unsubscribe, send an empty email to <a href="mailto:eeglablist-unsubscribe@sccn.ucsd.edu">eeglablist-unsubscribe@sccn.ucsd.edu</a><br>> For digest mode, send an email with the subject "set digest mime" to <a href="mailto:eeglablist-request@sccn.ucsd.edu">eeglablist-request@sccn.ucsd.edu</a><br>><br>><br>><br>> --<br>> Makoto Miyakoshi<br>> Swartz Center for Computational Neuroscience<br>> Institute for Neural Computation, University of California San Diego<br>> _______________________________________________<br>> Eeglablist page: <a href="http://sccn.ucsd.edu/eeglab/eeglabmail.html" target="_blank">http://sccn.ucsd.edu/eeglab/eeglabmail.html</a><br>> To unsubscribe, send an empty email to <a href="mailto:eeglablist-unsubscribe@sccn.ucsd.edu">eeglablist-unsubscribe@sccn.ucsd.edu</a><br>> For digest mode, send an email with the subject "set digest mime" to <a href="mailto:eeglablist-request@sccn.ucsd.edu">eeglablist-request@sccn.ucsd.edu</a><br><br></div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature"><div dir="ltr">Makoto Miyakoshi<br>Swartz Center for Computational Neuroscience<br>Institute for Neural Computation, University of California San Diego</div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Dec 18, 2014 at 4:32 PM, James Bonello <span dir="ltr"><<a href="mailto:jamesbonello9@gmail.com" target="_blank">jamesbonello9@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi EEGLablist,<div><br></div><div>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:</div><div><br></div><div><div>X = EEG.data(channel_no,:,epoch_no);</div><div> length_X = EEG.pnts; % signal length</div><div> nfft = 2^nextpow2(length_X); % Next power of 2 from length of x</div><div> </div><div> no_of_overlapped_samples = 0;</div><div> window = 512;</div><div> </div><div> % get psd using pwelch()</div><div> [spectra,freqs] = pwelch(X, window, no_of_overlapped_samples, nfft, EEG.srate);</div><div><br></div><div> spectra = 10*log(spectra);</div><div><br></div><div>Above I am calculating it using the pwelch() function, with spectopo it is as below:</div><div><br></div><div>[spectra,freqs] = spectopo(EEG.data(channel_no,:,epoch_no), 0, EEG.srate);<br></div><div><br></div><div><br></div><div>Is there something I am missing or am doing wrong?</div><div><br></div><div>Thank you!</div><span class="HOEnZb"><font color="#888888"><div><br></div>-- <br><div><b><i>James Bonello</i></b></div>
</font></span></div></div>
<br>_______________________________________________<br>
Eeglablist page: <a href="http://sccn.ucsd.edu/eeglab/eeglabmail.html" target="_blank">http://sccn.ucsd.edu/eeglab/eeglabmail.html</a><br>
To unsubscribe, send an empty email to <a href="mailto:eeglablist-unsubscribe@sccn.ucsd.edu">eeglablist-unsubscribe@sccn.ucsd.edu</a><br>
For digest mode, send an email with the subject "set digest mime" to <a href="mailto:eeglablist-request@sccn.ucsd.edu">eeglablist-request@sccn.ucsd.edu</a><br></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature"><div dir="ltr">Makoto Miyakoshi<br>Swartz Center for Computational Neuroscience<br>Institute for Neural Computation, University of California San Diego<br></div></div>
</div>