<div dir="ltr">Dear Samuel,<div><br></div><div>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.</div>
<div><br></div><div>-1 to 1 uV peak-to-peak sine wave measures 0dB at EEGLAB spectopo(), which makes sense to me.</div><div><br></div><div>Maktoo</div><div><br></div><div><br></div></div><div class="gmail_extra"><br><br>
<div class="gmail_quote">
On Wed, Jun 25, 2014 at 10:22 AM, Samuel Tomlinson <span dir="ltr"><<a href="mailto:stomlin1@swarthmore.edu" target="_blank">stomlin1@swarthmore.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div><div style="font-family:times new roman,new york,times,serif;font-size:12pt;color:#000000">Hello eeglablist,<br><div><br></div>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>
<div><br></div>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>
<div><br></div>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><div>
<br></div>Thank you,<br>Sam Tomlinson<br><div><br></div><div><br></div>Here's the code:<br><div><br></div>% 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><div><br></div>% get some params from EEG dataset<br>nchans = size(data,1);<br>N = size(data,2);<br>nepochs = size(data,3); <br><div><br></div>% set some standard values for fft<br>
nyquist = srate/2;<br>freqs = linspace(0,nyquist,floor(N/2+1));<br><div><br></div>% initialize output matrix<br>eegspec = zeros(nchans,floor(N/2+1));<br><div><br></div>% 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><div><br></div>% % 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><div><br></div>freqs = freqs(1:size(eegspec,2)); % make sure dimensions match<br><div><br></div><br>% ######## run comparison with EEGLAB spectopo ######<br><div><br></div>chan = randi([1,128]); % get a random channel<br>
<div><br></div>% 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><div><br></div>% first, show that frequencies are same<br>disp('frequency: ')<br>corrcoef(freqs,b)<br><div><br></div>% 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><div><br></div>%eegspec = eegspec.*sqrt(N);<br><div><br></div>% 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><div><br></div>% 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>
<div><br></div>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><div><br></div><br></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 dir="ltr">Makoto Miyakoshi<br>Swartz Center for Computational Neuroscience<br>Institute for Neural Computation, University of California San Diego<br></div>
</div>