tic
segmentDuration=52;
resampleRate=4096;
deltaF=0.25;
deltaT=1/resampleRate;
numPoints=segmentDuration*resampleRate;
psdFFTLength=resampleRate/deltaF;
psdWindow= hann(psdFFTLength,'periodic');
csdWindow= hann(numPoints,'periodic');
wind=csdWindow;
psdOverlapLength=psdFFTLength/2; 
detrendFlag='none';
%Choose the number of segments to analyze
numsegments=2000;
di=[];
dl=[];
%T1L=[];
%T2L=[];
%T1C=[];
%T2C=[];
%T1R=[];
%T2R=[];
%start doing what calCSD does
n1a=randn(numPoints/2,1);
n1b=randn(numPoints/2,1);
n1c=randn(numPoints/2,1);
n2a=randn(numPoints/2,1);
n2b=randn(numPoints/2,1);
n2c=randn(numPoints/2,1);
%[T1l,F]=psd([n1a;n1b],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag); 
%[T2l,F]=psd([n2a;n2b],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
%[T1c,F]=psd([n1b;n1c],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
%[T2c,F]=psd([n2b;n2c],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag); 
%This normalization makes fc as if from (0,1) data, and goes with sqrt(2*jsum) in snr.
nw=2/norm(wind)^2;
%This normalization accounts for frequency leakage, and goes with 1/sqrt(jsum) in snr.
%nw=1/numPoints/deltaF/2;
%Zero padding is done by windowAndFFT
zpad=zeros(numPoints,1);
%Prepare for coarseGrain.
numFreqs=1/2*resampleRate/deltaF;
cz=zeros(numFreqs,1);
jsum=2*deltaF*segmentDuration;
w=[1/2;ones(jsum-1,1);1/2];
outputsize=(numFreqs-1)*numsegments;
for in=1:numsegments
na=[wind.*[n1b;n1c];zpad];
nb=[wind.*[n2b;n2c];zpad];
fa=fft(na);
fb=fft(nb);
fa(1)=fa(1)+i*fa(end/2+1);
fb(1)=fb(1)+i*fb(end/2+1);
%This step is from calCSD.
fc=nw*conj(fa(1:end/2)).*fb(1:end/2);
%The normalization reflects the windowing and zero padding.
%Without them, the correct factor would be 4/numPoints.
%Start the coarseGrain calculation.
%determine statistics for independent points: std(real(d))=sqrt(52).
ci=cz;
for jn=1:jsum
ci=ci+fc((jn-1)*numFreqs+1:jn*numFreqs);
end
di=[di;ci];
%coarse grain adjacent points: notice the difference from above.
cl=cz(1:end-1);
for jn=1:jsum+1
cl=cl+w(jn)*fc(jsum/2+jn:jsum:end-3*jsum/2+jn);
end
dl=[dl;cl];
n1d=randn(numPoints/2,1);
n2d=randn(numPoints/2,1);
%[T1r,F]=psd([n1c;n1d],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
%[T2r,F]=psd([n2c;n2d],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
%T1L=[T1L,T1l];
%T2L=[T2L,T2l];
%T1C=[T1C,T1c];
%T2C=[T2C,T2c];
%T1R=[T1R,T1r];
%T2R=[T2R,T2r];
n1a=n1b;
n2a=n2b;
n1b=n1c;
n2b=n2c;
n1c=n1d;
n2c=n2d;
%T1l=T1c;
%T1c=T1r;
%T2l=T2c;
%T2c=T2r;
end
%T1L(1,:)=(T1L(1,:)+T1L(end,:))/2;
%T1C(1,:)=(T1C(1,:)+T1C(end,:))/2;
%T1R(1,:)=(T1R(1,:)+T1R(end,:))/2;
%T2L(1,:)=(T2L(1,:)+T2L(end,:))/2;
%T2C(1,:)=(T2C(1,:)+T2C(end,:))/2;
%T2R(1,:)=(T2R(1,:)+T2R(end,:))/2;
%Tal=reshape(T1L(2:end-1,:),outputsize,1);
%Tac=reshape(T1C(2:end-1,:),outputsize,1);
%Tar=reshape(T1R(2:end-1,:),outputsize,1);
%Ta=(Tal+Tar)/2;
%Tbl=reshape(T2L(2:end-1,:),outputsize,1);
%Tbc=reshape(T2C(2:end-1,:),outputsize,1);
%Tbr=reshape(T2R(2:end-1,:),outputsize,1);
%Tb=(Tbl+Tbr)/2;
%Tadj=sqrt(Ta.*Tb);
%Tloc=sqrt(Tac.*Tbc);
%Note: normalization here differs from coarseGrain output.
%dd=[real(d);imag(d)];
%snr=real(d)./Tadj/sqrt(jsum);
save('~/data/csddata.mat','di','dl');
%var(real(di))
%var(real(dl))
toc