cross-spectrum to cross-correlogram

[Jasper Poort]

Dear all,



I'd like to convert the cross-spectrum to a cross-correlogram. My first
approach would be simply to ifft the output of freqanalysis with output
option 'complex'.



However, I have difficulty understanding the result of this procedure and in
particular how I can get output comparable to the matlab function xcorr
using FieldTrip. Below I inserted a small test script.



Best, Jasper





% using the fake data from the FieldTrip tutorial 'Fourier analysis of

% oscillatory power and coherence'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%get many repetitions of two signals with somewhat consistent phase
difference

clear all

close all



frq = 10; % Hz

len = 1; % seconds

smpfrq = 100; % Hz

numrpt = 1000;

phsspreadfac = 0.5;

circulran = mod(randn(2.*numrpt,2).*phsspreadfac + pi, 2.* pi) - pi;

ranphs = circulran ./ (2 .* pi);

phsdif = 45 ./ 360;

noifac = 1./50;

for rptlop = 1:numrpt

    wav(:,rptlop,1) =
sin(((0:(len.*smpfrq-1))./(len.*smpfrq).*(frq.*2.*pi))+(ranphs(rptlop,1).*2.
*pi)) + ...

      randn(1,len.*smpfrq).*noifac;

    wav(:,rptlop,2) =
sin(((0:(len.*smpfrq-1))./(len.*smpfrq).*(frq.*2.*pi))+((ranphs(rptlop,2)+ph
sdif).*2.*pi)) + ...

      randn(1,len.*smpfrq).*noifac;

end



% get the FFT of the waves

for rptlop = 1:numrpt

    fftwav(:,rptlop,1) = fft(wav(:,rptlop,1));

    fftwav(:,rptlop,2) = fft(wav(:,rptlop,2));

end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% make FieldTrip structure

dum.label   = {'1','2'};

dum.fsample = smpfrq;

for i=1:size(wav,2), % for every trial

    dum.trial{i}=squeeze(wav(:,i,:))';

    dum.time{i} =[0:(len.*smpfrq-1)]./smpfrq;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%fft:

cfg.method     = 'mtmfft';

cfg.output     = 'fourier';

cfg.taper      = 'hanning';

cfg.tapsmofrq  = 1;

[freq] = freqanalysis(cfg, dum)



% cross-spectrum

CS=freq.fourierspctrm(:,1,:).*conj(freq.fourierspctrm(:,2,:));

y = real((squeeze(mean(CS))));

figure;plot(y)



% inverse fft of cross-spectrum

y = real(ifft(squeeze(mean(CS))));

figure;plot(y)



% compare to cross-correlation function

[c,lags] = xcorr(mean(squeeze(wav(:,:,1)),2),mean(squeeze(wav(:,:,2)),2));


figure;plot(lags,c)


----------------------------------
The aim of this list is to facilitate the discussion between users of the FieldTrip  toolbox, to share experiences and to discuss new ideas for MEG and EEG analysis. See also http://listserv.surfnet.nl/archives/fieldtrip.html and http://www.ru.nl/fcdonders/fieldtrip.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20080929/56b4d979/attachment-0001.html>