% eeglab2loreta
% needs the EEG. structure of EEGLAB
% downsamples to DOWN Hz
% calculates CSD cross-spectral density
% saves in LORETA-format
% 031114 sjo EEG. should only contain cephalic channels
% 030623 johannes.sarnthein@usz.ch with help of Roberto Pascual-Marqui
% 
DOWN=64; %downsample to DOWN Hz
WIN=2*DOWN; %2 seconds windowlength for continuous EEG
% WIN=EEG.pnts/EEG.srate*DOWN; % epochs of evoked EEG
NFFT=2^nextpow2(WIN);
clear csdmatrix
fprintf(1,'calculating multitaper cross spectral matrices .crs... \n')
for ic1=1:EEG.nbchan
    fprintf(1,' %i',ic1);
    for ic2=1:ic1
        data1=resample(double(EEG.data(ic1,:)),DOWN,EEG.srate);
        data2=resample(double(EEG.data(ic2,:)),DOWN,EEG.srate);
        [csdmatrix(ic2,ic1,:),freq]=csd(data1,data2,NFFT,DOWN,WIN,WIN/2); % so that (10,1)=0 and (1,10)=Cxy
    end 
end

header.NumFiles=1;
header.Ne=EEG.nbchan;
header.Nt=WIN;              % why is number of time frames needed for spectra ?
header.iFrq1=1;             % frequency limits are chosen in Loreta ...
header.iFrq2=length(freq);  % ... here the csdmatrix is described
header.NumMat=size(csdmatrix,3);
header.SampRateHz=DOWN;
header.frq1=freq(1);        % frequency limits are chosen in Loreta ...
header.frq2=freq(end);      % ... here the csdmatrix is described
header.FreqResHz=DOWN/WIN;
header.ClassicalBands='0';    % if true there are 7 or 8 csd matrices only
header.freq=freq;

fid=fopen('csd2loreta.crs','wb');

%write header
fwrite(fid,header.NumFiles,'int32');  % signed integer 32 bits = 4 byte
fwrite(fid,header.Ne,'int32');
fwrite(fid,header.Nt,'int32');
fwrite(fid,header.iFrq1,'int32');
fwrite(fid,header.iFrq2,'int32');
fwrite(fid,header.NumMat,'int32');
fwrite(fid,header.SampRateHz,'double');  % floating point 64 bits = 8 byte
fwrite(fid,header.frq1,'double');
fwrite(fid,header.frq2,'double');
fwrite(fid,header.FreqResHz,'double');      
fwrite(fid,header.ClassicalBands,'int8'); % one byte
%write all csdmatrices
for imat=1:header.NumMat % check at imat=12, i.e. 6 Hz
    realpart=real(csdmatrix(:,:,imat)); % leave real part in upper triangle
    imagpart=imag(csdmatrix(:,:,imat)).'; % transpose complex part to lower triangle
    matrix=realpart+imagpart;
    fwrite(fid,matrix,'double'); 
end
fclose(fid);
fprintf(1,'\n matrix written to file \n');

return